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O^'. Abstract 

2 ■ DDSCAT.6.0 is a freely available software package which applies the "discrete dipole approxima- 

I tion" (DDA) to calculate scattering and absorption of electromagnetic waves by targets with arbitrary 

■ geometries and complex refractive index. The DDA approximates the target by an array of polarizable 

points. DDSCAT.6.0 allows accurate calculations of electromagnetic scattering from targets with "size 
parameters" 27ra/A < 15 provided the refractive index m is not large compared to unity (|m — 1| < 1). 
K> ' DDSCAT.6.0 includes the option of using the FFTW (Fastest Fourier Transform in the West) package. 

^ . DDSCAT.6.0 also includes support for MP I (Message Passing Interface), permitting parallel calcula- 

tions on multiprocessor systems. 

The DDSCAT package is written in Fortran and is highly portable. The program supports calcu- 
lations for a variety of target geometries (e.g., ellipsoids, regular tetrahedra, rectangular solids, finite 
cylinders, hexagonal prisms, etc.). Target materials may be both inhomogeneous and anisotropic. It is 
straightforward for the user to "import" arbitrary target geometries into the code, and relatively straight- 
forward to add new target generation capabihty to the package. DDSCAT automatically calculates total 
cross sections for absorption and scattering and selected elements of the Mueller scattering intensity 
matrix for specified orientation of the target relative to the incident wave, and for specified scattering 
directions. 

This User Guide explains how to use DDSCAT.6.0 to carry out electromagnetic scattering calcula- 
tions. CPU and memory requirements are described. 
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1 Introduction 



DDSCAT.6.0 is a Fortran software package to calculate scattering and absorption of electromagnetic 
waves by targets with arbitrary geometries using the "discrete dipole approximation" (DDA). In this 
approximation the target is replaced by an array of point dipoles (or, more precisely, polarizable points); 
the electromagnetic scattering problem for an incident periodic wave interacting with this array of point 
dipoles is then solved essentially exactly. The DDA (sometimes referred to as the "coupled dipole 
approximation") was apparently first proposed by Purcell & Pennypacker (1973). DDA theory was 
reviewed and developed further by Draine (1988), Draine & Goodman (1993), and recently reviewed by 
Draine & Flatau (1994) and Draine (2000). 

DDSCAT.6.0 is a Fortran implementation of the DDA developed by the authors. It is intended to be 
a versatile tool, suitable for a wide variety of applications ranging from interstellar dust to atmospheric 
aerosols. As provided, DDSCAT.6.0 should be usable for many applications without modification, but 
the program is written in a modular form, so that modifications, if required, should be fairly straightfor- 
ward. 

The authors make this code openly available to others, in the hope that it will prove a useful tool. 
We ask only that: 

• If you publish results obtained using DDSCAT, please consider acknowledging the source of the 
code. 

• If you discover any errors in the code or documentation, please promptly communicate them to 
the authors. 

• You comply with the "copyleft" agreement (more formally, the GNU General Public License) of 
the Free Software Foundation: you may copy, distribute, and/or modify the software identified as 
coming under this agreement. If you distribute copies of this software, you must give the recip- 
ients all the rights which you have. See the file doc /copyleft distributed with the DDSCAT 
software. 

We also strongly encourage you to send email to the authors identifying yourself as a user of DDSCAT; 
this will enable the authors to notify you of any bugs, corrections, or improvements in DDSCAT. 

The current version, DDSCAT.6.0, uses the DDA formulae from Draine (1988), with dipole po- 
larizabilities determined from the Lattice Dispersion Relation (Draine & Goodman 1993). The code 
incorporates Fast Fourier Transform (FFT) methods (Goodman, Draine, & Flatau 1991). 

This User Guide assumes that you have already obtained the Fortran source code for DDSCAT.6.0 
either by download from |http : / /www, astro .princeton . edu/~draine |or by following the 
instructions in the README file.^ We refer you to the list of references at the end of this document for 
discussions of the theory and accuracy of the DDA [first see the recent reviews by Draine and Flatau 
(1994) and Draine (2000)]. In S|4]we describe the principal changes between DDSCAT.6.0 and the 

'To obtain the README file: 

(1) anonymous ftp to astro .princeton . edu , 

(2) cd draine/scat/DDA/ver5, and 

(3) get README. 
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previous releases.^ The succeeding sections contain instructions for: 

• compiling and linking the code; 

• running a sample calculation; 

• understanding the output from the sample calculation; 

• modifying the parameter file to do your desired calculations; 

• specifying target orientation; 

• changing the DIMENSIONing of the source code to accommodate your desired calculations. 

The instructions for compiling, linking, and running will be appropriate for a UNIX system; slight 
changes will be necessary for non-UNIX sites, but they are quite minor and should present no difficulty. 

Finally, the current version of this User Guide can be obtained from 
|http : //arxiv . orq/abs/astro-ph/0008151| - you will be offered the options of download- 
ing 

• Latex source 

• Postscript 

• Other formats - click on this to obtain the UserGuide as a PDF file. 



2 Applicability of the DDA 

The principal advantage of the DDA is that it is completely flexible regarding the geometry of the target, 
being limited only by the need to use an interdipole separation d small compared to ( 1 ) any structural 
lengths in the target, and (2) the wavelength A. Numerical studies (Draine & Goodman 1993; Draine & 
Flatau 1994; Draine 2000) indicate that the second criterion is adequately satisfied if 

\m\kd < 1 , (1) 

where m is the complex refractive index of the target material, and k = 2tt/X, where A is the wavelength 
in vacuo. However, if accurate calculations of the scattering phase function (e.g., radar or hdar cross 
sections) are desired, a more conservative criterion 

|m|fcd<0.5 (2) 

will ensure that differential scattering cross sections dCscn/dil are accurate to within a few percent of 
the average differential scattering cross section Csca/47r (see Draine 2000). 

Let V be the target volume. If the target is represented by an array of N dipoles, located on a cubic 
lattice with lattice spacing d, then 

V = Nd^ . (3) 
We characterize the size of the target by the "effective radius" 

aeff EE (3^^/4^)1/3 , (4) 

^The previous "official releases" were 

• DDSCAT.4b, 

• DDSCAT.4c -while never announced, DDSCAT . 4 c was made available to a number of interested users. 

• DDSCAT.5a7 

• DDSCAT.5a8 

• DDSCAT.5a9 

• DDSCAT.5alO 
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the radius of an equal volume sphere. A given scattering problem is then characterized by the dimen- 
sionless "size parameter" 

_ 27racff - 

X = fcacff = — — . (5) 



The size parameter can be related to N and \m\kd: 



27raoff 62.04 { N 
A 

Equivalently, the target size can be written 



1/3 



aeff = 9.873^ (^^j ■\m\kd . (7) 

Practical considerations of CPU speed and computer memory currently available on scientific worksta- 
tions typically limit the number of dipoles employed to < 10^ (see ill 61 for limitations on N due to 
available RAM); for a given N, the limitations on \m\kd translate into limitations on the ratio of target 
size to wavelength. 

For calculations of total cross sections Cabs and Csca, we require \ni\kd < 1: 

a^B < 9-88r^ { — ] oi x < { — ] ■ (8) 



106/ Iml VlO^ 



For scattering phase function calculations, we require \m\kd < 0.5: 

A f N 31.02 / A^ 

fleff < 4.94— — or X < —— T-TTT . (9) 



106 / Iml V 106 



It is therefore clear that the DDA is not suitable for very large values of the size parameter x, or very 
large values of the refractive index m. The primary utility of the DDA is for scattering by dielectric 
targets with sizes comparable to the wavelength. As discussed by Draine & Goodman (1993), Draine 
& Flatau (1994), and Draine (2000), total cross sections calculated with the DDA are accurate to a few 
percent provided A^ > 10^ dipoles are used, criterion Q is satisfied, and |m — 1| < 2. 

Examples illustrating the accuracy of the DDA are shown in Figs.[m2]which show overall scattering 
and absorption efficiencies as a function of wavelength for different discrete dipole approximations to a 
sphere, with N ranging from 304 to 59728. The DDA calculations assumed radiation incident along the 
(1,1,1) direction in the "target frame". Figs.|^}0show the scattering properties calculated with the DDA 
for X — ka — 7 . Additional examples can be found in Draine & Flatau (1994) and Draine (2000). 



3 DDSCAT.6.0 

3.1 What Does It Calculate? 

DDSCAT.6.0, like previous versions of DDSCAT, solves the problem of scattering and absorption by an 
array of polarizable point dipoles interacting with a monochromatic plane wave incident from infinity. 
DDSCAT.6.0 has the capability of automatically generating dipole array representations for a variety of 
target geometries (see ill9> and can also accept dipole array representations of targets supplied by the 
user (although the dipoles must be located on a cubic lattice). The incident plane wave can have arbitrary 
elliptical polarization (see ^2l\ . and the target can be arbitrarily oriented relative to the incident radiation 
(see ^n\ . The following quantities are calculated by DDSCAT.6.0 : 



6 




10 11 12 13 



Figure 1: Scattering and absorption for a sphere witli m = 1.33 + O.Oli. The upper panel shows the exact 
values of Qsca and Q^hs^ obtained with Mie theory, as functions of x = ka. The middle and lower panels 
show fractional errors in Qsca and Qahs, obtained using DDSCAT with polarizabilities obtained from the 
Lattice Dispersion Relation, and labelled by the number N of dipoles in each pseudosphere. After Fig. 1 of 
Draine & Flatau (1994). 
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Figure 2: Same as Fig.^ but for m = 2 + z. After Fig. 2 of Draine & Flatau (1994). 
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Figure 3: Differential scattering cross section for m = 1.33 + O.Oli pseudosphere and ka = 7. Lower 
panel shows fractional error compared to exact Mie theory result. The N = 17904 pseudosphere has 
\m\kd = 0.57, and an rms fractional error in da/dQ of 2.4%. After Fig. 5 of Draine & Flatau (1994). 




Figure 4: Same as Fig.|3lbut for m = 2 + i The = 59728 pseudosphere has \m\kd = 0.65, and an rms 
fractional error in da/dQ of 6.7%. After Fig. 8 of Draine & Flatau (1994). 
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• Absorption efficiency factor Qabs = Cabs/'i'aoff' where Cabs is the absorption cross section; 

• Scattering efficiency factor Qsca = Csca/Tra^g, where Csca is the scattering cross section; 

• Extinction efficiency factor Qcxt = Qsca + Qabsl 

• Phase lag efficiency factor Qpha, defined so that the phase-lag (in radians) of a plane wave after 
propagating a distance L is just ntQ-phsjia^gL, where nt is the number density of targets. 

• The 4x4 Mueller scattering intensity matrix Sij describing the complete scattering properties of 
the target for scattering directions specified by the user. 

• Radiation force efficiency vector Qj-ad (see i jl5> . 

• Radiation torque efficiency vector Qr (see JTst . 



3.2 Application to Targets in Dielectric Media 

Let ijj be the angular frequency of the incident radiation. For many applications, the target is essen- 
tially in vacuo, in which case the dielectric function e which the user should supply to DDSCAT is the 
actual complex dielectric function etargot('^), or complex refractive index mtargct(t^) = -Retarget of the 
target material. 

However, for many applications of interest (e.g., marine optics, or biological optics) the "target" 
body is embedded in a (nonabsorbing) dielectric medium, with (real) dielectric function emcdium(^), 
or (real) refractive index mmodium(t^) — \/emcdium- DDSCAT is fully apphcable to these scattering 
problems, except that: 

• The "dielectric function" or "refractive index" supplied to DDSCAT should be the relative dielec- 
tric function 

^medium \}^ } 

or relative refractive index: 

I ^ TOtargot(w) 

m[uj) = --. (11) 

Wnicdium(WJ 

• The wavelength A specified in ddscat . par should be the wavelength in the medium: 

A= , (12) 

^medium 

where Xyac — 27rc/a; is the wavelength in vacuo. 

The absorption, scattering, extinction, and phase lag efficiency factors Qabs, Qsca, and Qoxt calculated 
by DDSCAT will then be equal to the physical cross sections for absorption, scattering, and extinction 
divided by Tra^j^ - e.g., the attenuation coefficient for radiation propagating through a medium with 
a density rit of scatterers will be just a — ntQoxtTra^g . Similarly, the phase lag (in radians) after 



7ra 



2 

off- 



propagating a distance L will be n^Qphj 

The elements Sij of the 4x4 Mueller scattering matrix S calculated by DDSCAT will be correct for 
scattering in the medium: 



S • (13) 



\27rr J 

where Ii„ and Igca are the Stokes vectors for the incident and scattered Ught (in the medium), r is 
the distance from the target, and A is the wavelength in the medium (eg. I12>. See ^1231 for a detailed 
discussion of the Mueller scattering matrix. 

The time-averaged radiative force and torque (see ^\5\ on a target in a dielectric medium are 

Frad = QprTTOggWrad , (14) 
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2 ^ 

Trad = Qr7raeff'«rad2:^ , (15) 



where the time-averaged energy density is 



^rad — ^medium n ' (16) 

where Eq cos{u!t + (f>) is the electric field of the incident plane wave in the medium. 



4 What's New? 

DDSCAT.6.0 differs from DDSCAT.Sa in the following ways: 

1. Via the parameter file ddscat .par, the user can now select up to 9 elements of the MuUer 
scattering matrix Sij to be printed out. 

2. The maximum number of iterations allowed has been increased from 300 to 10000, since DDSCAT 
is now increasingly employed for computations that converge relatively slowly: large numbers of 
dipoles, large values of the scattering parameter, or large values of the refractive index. The num- 
ber of iterations used for a solution is recorded in the wwaarbbkccc . sea output files. 

3. A new target option is available: LYRSLB (a rectangular slab with up to 4 different material 
layers). 

4. User must now specify (in ddscat . par which IWAV, IRAD, IWAV to start with. The "default" 
choice will be 0, but when computations have been interupted by, e.g., power outage, there 
will be occasions when it is convenient to be able to resume with the first incomplete calculation. 

5. A new FFT option is supported: FFTW2 1. This invokes the FFTW ("Fastest Fourier Transform in 
the West") package of Frigo and Johnson. The calls are compatible with versions 2.1.x of FFTW. 
Instructions on compiling and linking to include FFTW are given below ( i|6.5t . 

6. FFT options BRENNR and TMPRTN have been eliminated as they seem to offer no advantages 
relative to GPFAFT or FFTW21. 

7. MP I is now supported. Users can now carry out parallel calculations using DDSCAT.6.0 on 
multiprocessor systems in conformance with the MP I (Message Passing Interface) standards ( ^24i . 
This of course requires that MP I be already installed on the system. 



5 Obtaining the Source Code 

The easiest way to download the source code is from the DDSCAT home page: 

|http : //www. astro .princeton . edu/^dralne/DDSCAT .html| 
where you can obtain the file ddscat 6 . . tar . gz - a "gzipped tarfile" containing the complete 
source code and documentation. This can also be obtained by anonymous ftp from astro, princeton. edu, 
directory draine/scat/ddscat/ver6 . 0. Note that it is a compressed (binary) file. 

The source code can be installed as follows. First, place the file ddscat 6.0. tar . gz in the 
directory where you would like DDSCAT to reside. You should have at least 5 Mbytes of disk space 
available. 

If you are on a Linux system, you should be able to type 
tar xvzf ddscat 6 . . tar . gz 
which will "gunzip" the tarfile and then "extract" the files from the tarfile. If your version of "tar" 
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doesn't support the "z" option (e.g., you are running Solaris) then try 

zcat ddscat 6 . . tar . gz | tar xvf - 
If neither of the above work on your system, try the two-stage procedure 

gunzip ddscat 6 . . tar . gz 

tar xvf ddscat 6 . . tar 
(The disadvantage of the two-stage procedure is that it uses more disk space, since after the second step 
you will have the uncompressed tarfile ddscat 6.0. tar - about 3.8 Mbytes - in addition to all the 
files you have extracted from the tarfile - another 4.6 Mbytes). 

Any of the above approaches should create subdirectories src, doc, misc, and IDL. The 
source code will be in subdirectory src, and documentation in subdirectory doc. 

If you are running Windows on a PC, you will need the "win zip" program, which can be down- 
loaded from http : / / www . win zip . com win zip should be able to "unzip" the gzipped tarfile 
dd scat 6.0. tar. gz and "extract" the various files from it automatically. 

6 Compiling and Linking 

In the discussion below, it will be assumed that the source code for DDSCAT.6.0 has been installed 
in a directory DDA/src. To compile the code on a Unix system, position yourself in the directory 
DDA/src. The Makefile supplied has compiler options set as appropriate for use of the g77 com- 
piler under the Linux operating system. If you have a different Fortran compiler, you will probably need 
to edit Makefile to provide the desired compiler options. If you are using an operating system other 
than Linux, you may need to change the Makefile choice of timeit. Makefile contains samples of 
compiler options for selected operating systems, including Linux, HP AUX, IBM AIX, and SGI IRIX 
operating systems. 

So far as we know, there are only two operating-system-dependent aspects of DDSCAT.6.0: (1) the 
device number to use for "standard output", and (2) the TIMEIT routine. There are, in addition, three 
installation-dependent aspects to the code: 

• the procedure for linking to the FFTW library (see ^6.5i 

• the procedure for linking to the netCDF library (see ill 0.31 for discussion of the possibility of 
writing binary files using the machine-independent net CDF standard). 

• the procedure for Unking to the MP I library (see ^24. 2\ . 

6.1 Device Numbers IDVOUT and IDVERR 

The variables IDVOUT and IDVERR specify device numbers for "running output" and "error mes- 
sages", respectively. Normally these would both be set to the device number for "standard output" (e.g., 
writing to the screen if running interactively). The variables IDVOUT and IDVERR are set by DATA 
statements in the "main" program DDSCAT. f and in the output routine WRIMSG (file wrimsg. f). 
Under Sun Fortran, DATA IDVOUT/ 0/ results in unbuffered output to "standard output"; unbuffered 
output (if available) is nice so that if you choose to direct your output to a file (e.g., using ddscat 
>& ddscat . out &) the output file will contain up-to-date information. Other operating systems or 
compilers may not support this, and you may need to edit DDSCAT . f to change the two data statements 
to DATA IDVOUT/ 6/ and DATA IDVERR/ 6 /, and edit wrimsg . f to change DATA IDVOUT/0/ 
to DATA IDVOUT/ 6/. 
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6.2 Subroutine T I ME I T 

The only other operating system-dependent part of DDSCAT.6.0 is the single subroutine TIMEIT. 
Several versions of TIMEIT are provided: 

• timeit_sun . f uses the SunOS system call etime 

• timeit-convex . f uses the Convex OS system call etime 

• timeit_cray . f uses the system call second 

• timeit_hp . f uses the HP-AUX system calls sysconf and times 

• timeit_ibm6000 . f uses the AIX system call mclock 

• timeit_osf . f uses the DEC OSF system call etime 

• timeit_sgi . f uses the IRIX system call etime 

• timeit_vms . f uses the VMS system calls LIB$ INIT_TIMER and LIB$SHOW_TIMER 

• timeit_titan . f uses the system call cputim 

• timeit_null . f is a dummy routine which provides no timing information, but can be used 
under any operating system. 

You must compile and Hnk one of the timeitjacx . f routines, possibly after modifying it to work on 
your system; timeit_null . f is the easiest choice, but it will return no timing information. 

6.3 Optimization 

The performance of DDSCAT.6.0 will benefit from optimization during compilation and the user should 
enable the appropriate compiler flags (e.g., g77 -c -0, or pgf77 -c -fast). 

There is one exception: the LAPACK routine SLAMCl, which is called to determine the number of 
bits of precision provided by the computer architectures, should NOT be optimized, because under some 
optimizers the resulting code will end up in an endless loop. Therefore the Makefile has a separate com- 
pilation flag fFLAGSslamcl which should not specify optimization (e.g., just g7 7 -c, or pgf 7 7 
-c).4 

6.4 Leaving FFTW Capability Disabled 

FFTW^ - "Fastest Fourier Transform in the West" - is a publicly available FFT implementation that 
performs very well (see ^13\ . DDSCAT.6.0 supports use of FFTW, but it is recommended that first-time 
users of DDSCAT first get the code up and running without FFTW support, using the built-in GPFA FFT 
routine written by Clive Temperton, which performs almost as well as FFTW. This is also the option 
you will have to follow if FFTW is not installed on your system (though you can install it yourself, if 
necessary - see ii6.5l below). 

The Makefile supplied with the DDSCAT.6.0 distribution is initially set up to leave support for FFTW 
disabled by using a "dummy" subroutine dummycxf f tw . f . 

When you type make ddscat you will create a version of ddscat which will not recognize 
option FFTW2 1 - you must specify option GFPAFT. 

^Non-Unix sites: The VMS-compatible version of TIMEIT is included in the file SRC 9 . FOR. For non-VMS sites, you will 
need to replace this version of TIMEIT with one of the other versions, which are to be found in the file MI SC. FOR. If you 
are having trouble getting this to work, choose the "dummy" version of TIMEIT from MISC .FOR: this will return no timing 
information, but is platform-independent. 

''Note that the LAPACK routines are only used for some target geometries, and even when used the computation time is insignif- 
icant. 

" Ihttp : / / www . f f tw. orgi 
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6.5 Enabling FFTW Capability 



The FFTW ("Fastest Fourier Transform in the West") package of Matteo Frigo and Steven G. Johnson 
appears to be one of the fastest public-domain FFT packages available (see ill3l for a performance 
comparison), and DDSCAT.6.0 allows this package to be used. In order to support the FFTW2 1 option, 
however, it is necessary to first install the FFTW package on your system (FFTW 2.1.x - we have not 
yet implemented support for FFTW 3.0, which has a different interface). 
FFTW offers two advantages: 

• It appears, for some cases, to be somewhat faster than the GPFA FFT algorithm which is akeady 
included in the DDSCAT package (see Jojand Fig.|5}. 

• It does not restrict the dimensions of the "computational volume" to be factorizable as 2'3-'5'^, 
and therefore for some targets will use a smaller computational volume (and therefore require less 
memory). 

If the FFTW package is not yet installed on your system, then 

1. Download the FFTW 2.1.x package from |http : / /www . f f tw . orq| (as of this writing the 
latest version 2.1.x is f f tw-2 . 1 . 5 . tar . gz) into some convenient directory. 

2. tar xvfz f f tw-2 . 1 . 5 . tar . gz 

3. cd fftw-2.1.5 

4. ./configure — enable-float — enable-type-pref ix — prefix=PATH 
where PATH is the path to the directory in which you would like the fftw library installed (you 
must have write permission, of course). If you choose to omit 

— pref ix=PATH 

then fftw will be installed in the default directory /usr/local/, but unless you are root you 
may not have write permission to this directory. 

5. make 

6. make install 

This should install the file libsf f tw . a in the directory /usr/local/lib (or PATH/lib if you 
used option — pref ix=PATH when running configure). 
Once FFTW is installed, you will need to edit the Makefile: 

• define cxf ftw = cxf ftw 

• define LIBFFTW to give the correct path to libsf ftw . a 

Once the above is done, typing make ddscat should create an executable ddscat which will 
recognize option FFTW21 in the file ddscat .par. Note that you will be Unking fortran to C, and 
different compilers may behave idiosyncratically. The Makefile has examples which work for g 7 7 and 
pgf 7 7. If you have trouble, consult with someone familiar with Unking fortran and C modules on your 
system. 

6.6 Leaving the netCDF Capability Disabled 

netCDF^ is a standard for data transport used at many sites (see i ll0.3l for more information on net CDF). 
The Makefile supplied with the distribution of DDSCAT.6.0 is set up to link to a "dummy" subrou- 
tine dummywritenet . f instead of subroutine writenet . f , in order to minimize problems during 
initial compilation and linking. The "dummy" routine has no functionality, other than bailing out with a 
fatal error message if the user makes the mistake of trying to specify one of the netCDF options (ALLCDF 

^http://www.unidata.ucasr.edu/packages/netcdf/ 
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or ORICDF). First-time users of DDSCAT.6.0 should not try to use the netCDF option - simply skip 
this section, and specify option NOTCDF in ddscat . par. After successfully using DDSCAT.6.0 you 
can return to il6.7l to try to enable the netCDF capability. 

6.7 Enabling the netCDF Capability 

Subroutine WRITENET (file writenet . f) provides the capability to output binary data in the netCDF 
standard format (see iilQ.3t . In order to use this routine (instead of dummywritenet . f), it is neces- 
sary to take the following steps;^ 

1. Have netCDF already installed on your system (check with your system administrator). 

2. Find out where netcdf . inc is located and edit the include statement in writenet . f to 
specify the correct path to netcdf. inc. 

3. Find out where the libnetcdf.a library is located, and edit the Makefile so that the variable 
LIBNETCDF specifies the coiTect path to this library. 

4. Edit the Makefile to "comment out" (with a # symbol in column 1) the line writenet = 
dummywritenet and "uncomment" (remove the # symbol) the line writenet = writenet 
so that writenet . f will be compiled instead of the dummy routine dummywritenet . f . 

6.8 Compiling and Linking... 

After suitably editing the Makefile (while still positioned in DDA/src) simply type** 

make ddscat 

which should create an executable file DDA/ src/ ddscat . If the default (as-provided) Makefile is 
used, the g7 7 fortran compiler will be used to create an executable which: 

• will not have MP I support; 

• will nof have FFTW support; 

• will not have netCDF capability; 

• will contain timing instructions compatible with Linux, Solaris 1.x or 2.x, as well as several other 
versions of Unix. 

To add FFTW capability, see %.5\ To add netCDF capability, see ^6.7\ To add MP I support, see 
^241 To replace the Linux-compatible and Sun-compatible timing routine with another, see il6.2l 

6.9 Installation on Non-Unix Systems 

DDSCAT.Sa is written in standard Fortran-77 plus the DO ... ENDDO extension which appears to 
be supported by all current Fortran-77 compatible compilers. It is possible to run DDSCAT on non-Unix 
systems. If the Unix "make" utility is not available, here in brief is what needs to be accomplished: 

All of the necessary Fortran code to compile and link DDSCAT.Sa is included in the follow- 
ing files: SRCO . FOR, SRCl . FOR, SRC2 . FOR, SRC3 . FOR, SRC4 . FOR, SRC5 . FOR, SRC6 . FOR, 
SRC7 .FOR, SRC8 .FOR, SRC9 .FOR, CGCOMMON . FOR, GPFA.FOR, LAPACK.FOR, andPIM.FOR. 
There is an additional file MISC . FOR, but this is not needed for "basic" use of the code (see below). 

'Non-UNIX sites: You need to replace the "dummy" version of SUBROUTINE WRITENET in SRCl . FOR with the version 
provided in MISC . FOR. You will also need to consult your systems administrator to verify that the netCDF library has already 
been installed on your system, and to find out how to link to the library routines. 

•^Non-UNIX sites: see 
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The main program DDSCAT is found in SRCO .FOR. It calls a number of subroutines, which are 
included in the other * . FOR files. 

Four of the subroutines exist in more than one version. The "default" version of each is located in 
the file SRC 1 . FOR: 

• Selectthe version of SUBROUTINE CXFFTWfrom SRCl . FOR instead of the version in MI SC . FOR. 
The resulting code will not support FFTW capability. If you later wish to add FFTW capability, 
you will first need to install the FFTW library on your system (see ^6.5> . After you have done so, 
you can use the version of SUBROUTINE CXFFTW provided in MISC . FOR. 

• Select the version of SUBROUTINE WRITENET from SRCl .FOR instead of the version in 
MI SC. FOR. The resulting code will not support netCDF capability. (If netCDF capability is 
required, you will need to install netCDF libraries on your system - see %.7i . 

• Selectthe version of SUBROUTINE C3DFFT from SRCl . FOR (do not use the version in MISC . FOR) 

• There is one version of SUBROUTINE TIMEIT included in SRCl . FOR, and a number of addi- 
tional versions in MISC . FOR. Use the version from SRCl . FOR, which begins 

SUBROUTINE TIME IT (CMSGTM, DTIME ) 

C 

C timeit_null 
C 

C This version of timeit is a dummy which does not provide any 
C timing information. 

This version does not use any system calls, and therefore the code should compile and link without 
problems; however, when you run the code, it will not report any timing information reporting how 
much time was spent on different parts of the calculation. 

If you wish to obtain timing information, you will need to find out what system calls are supported 
by your operating system. You can look at the other versions of SUBROUTINE TIMEIT in 
MISC. FOR to see how this has been done under VMS and various version of Unix. 

Once you have selected the appropriate versions of CXFFTW, WRITENET, C3DFFT, and TIMEIT, you 
can simply compile and link just as you would with any other Fortran code with a number of modules. 
You should end up with an executable with a (system-dependent) name like DDSCAT .EXE. 

In addition to program DDSCAT, we provide two other Fortran 77 programs which may be useful. 
Program CALLTARGET can be used to call the target generation routines which may be helpful for 
testing purposes. Program TSTFFT is useful for comparing speeds of different FFT options. 

7 Moving the Executable 

Now reposition yourself into the directory DDA (e.g., type cd . .), and copy the executable from 
src/ddscat to the DDA directory by typing 

cp src/ddscat ddscat 
This should copy the file DDA/ src/ ddscat to DDA/ ddscat (you could, of course, simply create a 
symbolic link instead). Similarly, copy the sample parameter file ddscat . par and the file diel . tab 
to the DDA directory by typing 

cp doc/ddscat . par ddscat. par 

cp doc/diel.tab diel. tab 
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8 The Parameter File ddscat.par 



The directory DDA should now contain a sample file ddscat .par which provides parameters to the 
program ddscat. As provided (see Appendi?jXt. the file ddscat . par is set up to calculate scattering 
by a 8x6x4 rectangular array of 192 dipoles, with an effective radius — Ifim, at a wavelength of 
6.2832/im (for a "size parameter" 27raoft/A = 1). 

The dielectric function of the target material is provided in the file diel . tab, which is a sample 
file in which the refractive index is set to m = 1.33 + O.Olz at all wavelengths; the name of this file is 
provided to ddscat by the parameter file ddscat . par. 

The sample parameter file as supplied calls for the GPFA FFT routine (GPFAFT) of Temperton 
(1992) to be employed and the PBCGST iterative method to be used for solving the system of linear 
equations. (See section ill 31 and ill 41 for discussion of choice of FFT algorithm and choice of equation- 
solving algorithm.) 

The sample parameter file specifies (via option LATTDR) that the "Lattice Dispersion Relation" 
(Draine & Goodman 1993) be employed to determine the dipole polarizabilities. See ill ll for discussion 
of other options. 

The sample parameter files specifies options NOTE IN and NOTCDF so that no binary data files and 
no NetCDF files will be created by ddscat. 

The sample ddscat . par file specifies that the calculations be done for a single wavelength (6.2832/im) 
and a single effective radius (1.00/im). Note that in DDSCAT.6.0 the "effective radius" Ooff is the ra- 
dius of a sphere of equal volume - i.e., a sphere of volume Ncfi , where d is the lattice spacing and 
N is the number of occupied (i.e., non-vacuum) lattice sites in the target. Thus the effective radius 
Oeff = (37V/47r)i/3d . 

The incident radiation is always assumed to propagate along the x axis in the "Lab Frame". The sam- 
ple ddscat . par file specifies incident polarization state Gqi to be along the y axis (and consequently 
polarization state eo2 will automatically be taken to be along the z axis). I0RTH=2 in ddscat . par 
calls for calculations to be carried out for both incident polarization states (gqi and eo2 - see i|21> . 

The target is assumed to have two vectors ai and sl2 embedded in it; a2 is perpendicular to ai . In 
the case of the 8x6x4 rectangular array of the sample calculation, the vector ai is along the "long" 
axis of the target, and the vector a.2 is along the "intermediate" axis. The target orientation in the Lab 
Frame is set by three angles: /3, Q, and $, defined and discussed below in i|17l Briefly, the polar angles 

8 and $ specify the direction of ai in the Lab Frame. The target is assumed to be rotated around ai 
by an angle p. The sample ddscat . par file specifies /? = and $ = (see lines in ddscat . par 
specifying variables BETA and PHI), and calls for three values of the angle 8 (see line in ddscat . par 
specifying variable THETA). DDSCAT.6.0 chooses 8 values uniformly spaced in cos 8; thus, asking 
for three values of 8 between and 90° yields 8 = 0, 60°, and 90°. 

AppendixlAlprovides a detailed description of the file ddscat . par.^ 

9 Running DDSCAT.6.0 Using the Sample dds cat . par File 

To execute the program on a UNIX system (running either sh or csh), simply type 

ddscat >& ddscat. out & 
which will redirect the "standard output" to the file ddscat . out, and run the calculation in the back- 
ground. The sample calculation (8x6x4=192 dipole target, 3 orientations, two incident polarizations, 
with scattering calculated for 14 distinct scattering directions), requires 1.0 cpu seconds on a Sun Ultra- 
170 workstation. 

'Note that the structure of dds cat .par depends on the version of DDSCAT being used. Make sure you update old parameter 
tiles before using them with DDSCAT.6.0 ! 
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10 Output Files 



10.1 ASCII files 

If you run DDSCAT using the command 
ddscat >& ddscat.out & 
you will have various types of ASCII files when the computation is complete: 

• a file ddscat . out; 

• afilemtable; 

• a file qtable; 

• a file qtable2; 

• files wxxryyori . avg (one, wOOrOOori. avg, for the sample calculation); 

• if ddscat . par specified IWRKSC=1, there will also be files vxxryy'kzzz . sea (3 for the sample 
calculation: wOOrOOkOOO . sea, wOOrOOkOOl . sea, w00r00k002 . sea). 

The file ddscat . out will contain any error messages generated as well as a running report on the 
progress of the calculation, including creation of the target dipole array. During the iterative calcula- 
tions, Qoxt, Qabs, and Qpha are printed after each iteration; you will be able to judge the degree to 
which convergence has been achieved. Unless TIME IT has been disabled, there will also be timing 
information. 

The file mtable contains a summary of the dielectric constant used in the calculations. 

The file qtable contains a summary of the orientationally-averaged values of Qcxt, Qabs, Qsca, 
g{l) — {cos{6s)), and Q^k- Here Qoxt, Qabs, and Qsca are the extinction, absorption, and scattering 
cross sections divided by Tra^g. Qbk is the differential cross section for backscattering (area per sr) 
divided by Tra^jj . 

The file qtable2 contains a summary of the orientationally-averaged values of Qpha, QpoU and 
Qcpoi- Here Qpha is the "phase shift" cross section divided by Tra^fj (see definition in Draine 1988). 
Qpoi is the "polarization efficiency factor", equal to the difference between Qcxt for the two orthogonal 
polarization states. We define a "circular polarization efficiency factor" Qcpoi = QpoiQpha, since an 
optically-thin medium with a small twist in the alignment direction will produce circular polarization in 
initially unpolarized light in proportion to Qcpoi- 

For each wavelength and size, DDSCAT.6.0 produces a file with a name of the form 
vixxryyori . avg, where index xx (=00, 01, 02....) designates the wavelength and index yy (=00, 01, 
02...) designates the "radius"; this file contains Q values and scattering information averaged over how- 
ever many target orientations have been specified (see The file wOOrOOori.avg produced by the 
sample calculation is provided below in AppendixiBl 

In addition, if ddscat . par has specified IWRKSC=1 (as for the sample calculation), DDSCAT.6.0 
will generate files with names of the form wxxryykzzz . avg, where xx and yy are as before, and index 
zzz =(000,001,002...) designates the target orientation; these files contain Q values and scattering infor- 
mation for each of the target orientations. The structure of each of these files is very similar to that of 
the vixxryyori . avg files. Because these files may not be of particular interest, and take up disk space, 
you may choose to set IWRKSC=0 in future work. However, it is suggested that you run the sample 
calculation with IWRKSC=1. 

The sample ddscat . par file specifies IWRKSC=1 and calls for use of 1 wavelength, 1 target size, 
and averaging over 3 target orientations. Running DDSCAT.6.0 with the sample ddscat .par file 
will therefore generate files wOOrOOkOOO . sea, wOOrOOkOOl . sea, and w00r00k002 . sea . To 
understand the information contained in one of these files, please consult Appendix Icl which contains 
an example of the file wOOrOOkOOO. sea produced in the sample calculation. 
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10.2 Binary Option 



It is possible to output an "unformatted" or "binary" file (dd.bin) with fairly complete information, 
including header and data sections. This is accomplished by specifying either ALLBIN or ORIBIN in 
ddscat . par . 

Subroutine writebin. f provides an example of how this can be done. The "header" section 
contains dimensioning and other variables which do not change with wavelength, particle geometry, 
and target orientation. The header section contains data defining the particle shape, wavelengths, parti- 
cle sizes, and target orientations. If ALLBIN has been specified, the "data" section contains, for each 
orientation, Mueller matrix results for each scattering direction. The data output is limited to actual di- 
mensions of arrays; e.g. ns cat , 4 , 4 elements of Muller matrix are written rather than mxs cat , 4 , 4 . 
This is an important consideration when writing postprocessing codes. 

A skeletal example of a postprocessing code was written by us (PJF) and is provided in subdirectory 
DDA/IDL. If you do plan to use the Interactive Data Language (IDL) for postprocessing, you may 
consider using the netCDF binary file option which offers substantial advantages over the FORTRAN 
unformatted write. More information about IDL is available at http: / /www . rsinc . com/idl. 
Unfortunately IDL requires a license and hence is not distributed with DDSCAT. 

10.3 Machine-Independent Binary File Option: netCDF 

The "unformatted" binary file is compact, fairly complete, and (in comparison to ASCII output files) is 
efficiently written from FORTRAN. However, binary files are not compatible between different com- 
puter architectures if written by "unformatted" WRITE from FORTRAN. Thus, you have to process the 
data on the same computer architecture that was used for the DDSCAT calculations. We provide an 
option of using netCDF with DDSCAT. The netCDF library'" defines a machine-independent format for 
representing scientific data. Together, the interface, library, and format sup port the creation, access, and 
sharing of scientific data. For more information, go to the netcdf websitai^. 

Several major graphics packages (for example IDL) have adopted netCDF as a standard for data 
transport. For this reason, and because of strong and free support of netCDF over the network by 
UNIDATA, we have implemented a netCDF interface in DDSCAT. This may become the method of 
choice for binary file storage of output from DDSCAT. 

After the initial "learning curve" netCDF offers many advantages: 

• It is easy to examine the contents of the file (using tools such as ncdump). 

• Binary files are portable - they can be created on a supercomputer and processed on a workstation. 

• Major graphics packages now offer netCDF interfaces. 

• Data input and output are an order of magnitude faster than for ASCII files. 

• Output data files are compact. 
The disadvantages include: 

• Need to have netCDF installed on your system. 

• Lack of support of complex numbers. 

• Nontrivial learning curve for using netCDF. 

The public-domain netCDF software is available for many operating systems from 
|http : /'/www . unid ata . uca r . edu/packaqe s/netcdf . The steps necessary for enabling the 
netCDF capability in DDSCAT.6.0 are listed above in i|6.7l Once these have been successfully ac- 
complished, and you are ready to run DDSCAT to produce netCDF output, you must choose either 

' |http : / /www . unidata . ucar . edu /packages/ net cdf/| 
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the ALLCDF or ORICDF option in ddscat . par; ALLCDF will result in a file which will include the 
MuUer matrix for every wavelength, particle size, and orientation, whereas ORICDF will result in a file 
hmited to just the orientational averages for each wavelength and target size. 

11 Dipole Polarizabilities 

Option LATTDR specifies that the "Lattice Dispersion Relation" (Draine & Goodman 1993) be em- 
ployed to determine the dipole polarizabilities.; at this time this is the only supported option. 

12 Dielectric Functions 

In order to assign the appropriate dipole polarizabilities, DDSCAT.6.0 must be given the dielectric 
constant of the material (or materials) of which the target of interest is composed. Unless the user wishes 
to use the dielectric function of either solid or liquid H2O (see below), his information is supplied to 
DDSCAT through a table (or tables), read by subroutine DIELEC in file dielec . f, and providing 
either the complex refractive index m — n + ik 01 complex dielectric function e = ei + ie2 as a function 
of wavelength A. Since m — e^/^, or e = m^, the user must supply either m or e, but not both. 

The table formatting is intended to be quite flexible. The first line of the table consists of text, up 
to 80 characters of which will be read and included in the output to identify the choice of dielectric 
function. (For the sample problem, it consists of simply the statement m = 1.33 + O.Oli.) The 
second line consists of 5 integers; either the second and third or the fourth and fifth should be zero. 

• The first integer specifies which column the wavelength is stored in. 

• The second integer specifies which column Re(m) is stored in. 

• The third integer specifies which column lm(r7i) is stored in. 

• The fourth integer specifies which column Re(e) is stored in. 

• The fifth integer specifies which column lm(e) is stored in. 

If the second and third integers are zeros, then DIELEC will read Re(e) and lm(e) from the file; if the 
fourth and fifth integers are zeros, then Re(TO) and lm(m) will be read from the file. 

The third line of the file is used for column headers, and the data begins in line 4. There must be at 
least 3 lines of data: even if e is required at only one wavelength, please supply two additional "dummy" 
wavelength entries in the table so that the interpolation apparatus will not be confused. 

In the event that the user wishes to compute scattering by targets with the refractive index of either 
solid or liquid H2O, we have included two "built-in" dielectric functions. If H20ICE is specified on line 
10 of ddscat . par, DDSCAT will use the dielectric function of H2O ice at T=250K as compiled by 
Warren (1984). If H20LIQ is specified on line 10 of ddscat .par, DDSCAT will use the dielectric 
function for liquid H2O at T=280K using subroutine REFWAT written by Eric A. Smith. For more 
information see http://atol.ucsd. edu/^pf latau/ref rtab/index . htm| . 

13 Choice of FFT Algorithm 

One major change in going from DDSCAT .4b to 4c was modification of the code to permit use of the 
GPFA FFT algorithm developed by Dr. Clive Temperton (1992). In going from DDSCAT.SalO to 6.0 
we introduce a new FFT option: the FFTW package of Frigo and Johnson. DDSCAT continues to offer 
Temperton's GPFA code as an alternative FFT implementations. 
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Figure 5: Comparison of cpu time required by 3 different FFT implementations. It is seen that the GPFA 
and FFTW implementations have comparable speeds, much faster than Brenner's FFT implementation. 



To compare the performance of the GPFA algorithm and the FFTW code, we provide a program 
TSTFFT to compare the performance of different 3-D FFT implementations. The TSTFFT code re- 
quires the FFTW library (see %.5\ . TSTFFT also provides timings for an FFT implementation due to 
Brenner (included in previous releases of DDSCAT). 

To compile, link, and run this program on a Unix system," position yourself in the DDA/src 
directory and type 

make tstfft 
tstf ft 

Output results will be written into a file res . dat. TableQ]is the res . dat file created when run on a 
2.0 GHz Intel Xeon system, using the pgf77 compiler: 

It is clear that both the GPFA code and the FFTW code are much faster than the Brenner FFT (by 
a factor of 7 for the lOOx lOOx 100 case). The FFTW code and GPFA code are quite comparable in 
performance - for some cases the GPFA code is faster, for other cases the FFTW code is faster For 
target dimensions which are factorizable as 2*3'' 5*^ (for integer i, j, k), the GPFA and FFTW codes 
have the same memory requirements. For targets with extents Nx, Ny, Nz which are not factorizable 
as 2*3^5*^, the GPFA code needs to "extend" the computational volume to have values of Nx, Ny, and 
Nz which are factorizable by 2, 3, and 5. For these cases, GPFA requires somewhat more memory than 
FFTW. However, the fractional difference in required memory is not large, since integers factorizable as 
2*3-' 5*^ occur fairly frequently.'^ 

"Non-Unix systems: the TSTFFT Fortran source code is in the file MISC . FOR. 

''2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36, 40, 45, 48, 50, 54, 60, 64, 72, 75, 80, 81, 90, 96, 100, 108, 120, 
125, 128, 135, 144, 150, 160, 162, 180, 192, 200 are the integers < 200 whicli are of the form 2*3^ 5*=. 
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Table 1: FFT timing output from program tstf f t 



CPU time (sec) for different 3-D FFT methods 
Machine=2 . GHz Xeon, 256kB L2 cache, pgf77 compiler 
parameter LVR = 64 

LVR="length of vector registers" in gpf a2f , gpf a3f , gpf a5f 
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Based on Figure 5, we see that while for some cases FFTW 2 . 1 . 5 is faster than the GPFA algo- 
rithm, the difference appears to be fairly marginal. 

The GPFA code contains a parameter LVR which is set in data statements in the routines gpf a2 f , 
gpf a3f, and gpf a5f . LVR is supposed to be optimized to correspond to the "length of a vector reg- 
ister" on vector machines. As delivered, this parameter is set to 64, which is supposed to be appropriate 
for Grays other than the C90 (for the G90, 128 is supposed to be preferable). The value of LVR is not 
critical for scalar machines, as long as it is fairly large. We found little difference between LVR=64 and 
128 on a Sparc 10/51, on an Ultrasparc 170, and on an Intel Xeon cpu. You may wish to experiment 
with different LVR values on your computer architecture. To change LVR, you need to edit gpf a . f 
and change the three data statements where LVR is set. 

The choice of FFT implementation is obtained by specifying one of: 

• FFTW21 to use the FFTW 2.1.x algorithm of Frigo and Johnson - this is recommended, but 
requires that the FFTW 2.1.x library be installed on your system. 

• GPF AFT to use the GPFA algorithm (Temperton 1992) - this is almost as fast as the FFTW 
package, and does not require that the FFTW package be installed on your system; 

• CONVEX to use the native FFT routine on a Convex system (this option is presumably obsolete, 
but is retained as an example of how to use a system-specific native routine). 



14 Choice of Iterative Algorithm 

As discussed elsewhere (e.g., Draine 1988), the problem of electromagnetic scattering of an incident 
wave Einc by an array of point dipoles can be cast in the form 

AP = E (17) 

where E is a 3A^-dimensional (complex) vector of the incident electric field Ejnc at the N lattice sites, 
P is a 37V-dimensional (complex) vector of the (unknown) dipole polarizations, and A is a 3A^ x 3A^ 
complex matrix. 

Because 3A^ is a large number, direct methods for solving this system of equations for the unknown 
vector P are impractical, but iterative methods are useful: we begin with a guess (typically, P — 0) for 
the unknown polarization vector, and then iteratively improve the estimate for P until equation ( I17> is 
solved to some error criterion. The error tolerance may be specified as 

lAtAP-AtEl , 

- <h , (18) 



|AtE| 

where A^ is the Hermitian conjugate of A [(ylt)y = (Aji)*], and h is the error tolerance. We typically 
use h = 10~^ in order to satisfy eaAlli to high accuracy. The error tolerance h can be specified by the 
user (see AppendixIXl. 

A major change in going from DDSCAT.4b to 5a (and subsequent versions) was the implemen- 
tation of several different algorithms for iterative solution of the system of complex linear equations. 
DDSCAT.Sa and DDSCAT.6.0 have been structured to permit solution algorithms to be treated in a 
fairly "modular" fashion, facilitating the testing of different algorithms. Many algorithms were com- 
pared by Flatau (1997)''*; two of them performed well and are made available to the user in DDSCAT.6.0. 
The choice of algorithm is made by specifying one of the options: 



'^^In place of "preferable" users are encouraged to read "necessary" - there is some basis for fearing that results computed on a 
C90 with LVR other than 128 run the risk of being incorrect ! Please use LVR= 1 2 8 if running on a C90 ! 
'*A postscript copy of this report - file eg . ps - is distributed with the DDSCAT.6.0 documentation. 
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• PBCGST - Preconditioned BiConjugate Gradient with STabilitization method from the Parallel 
Iterative Methods (PIM) package created by R. Dias da Cunha and T. Hopkins. 

• PETRKP - the complex conjugate gradient algorithm of Petravic & Kuo-Petravic (1979), as coded 
in the Complex Conjugate Gradient package (CCGPACK) created by P.J. Flatau. This is the 
algorithm discussed by Draine (1988) and used in previous versions of DDSCAT. 

Both methods work well. Our experience suggests that PBCGST is often faster than PETRKP, by perhaps 
a factor of two. We therefore recommend it, but include the PETRKP method as an alternative. 

The Parallel Iterative Methods (PIM) by Rudnei Dias da Cunha (rdd@ukc . ac . uk) and Tim Hop- 
kins (trh@ukc.ac.uk) is a collection of Fortran 77 routines designed to solve systems of linear 
equations on parallel and scalar computers using a variety of iterative methods (available at 
|http: / /www, mat .ufrqs .br/pim-e .html^ . PIM offers a number of iterative methods, includ- 
ing 

• Conjugate-Gradients, CG (Hestenes 1952), 

• Bi-Conjugate-Gradients, BICG (Fletcher 1976), 

• Conjugate-Gradients squared, CGS (Sonneveld 1989), 

• the stabilised version of Bi-Conjugate-Gradients, BICGSTAB (van der Vorst 1991), 

• the restarted version of BICGSTAB, RBICGSTAB (Sleijpen & Fokkema 1992) 

• the restarted generalized minimal residual, RGMRES (Saad 1986), 

• the restarted generalized conjugate residual, RGCR (Eisenstat 1983), 

• the normal equation solvers, CGNR (Hestenes 1952 and CGNE (Craig 1955), 

• the quasi-minimal residual, QMR (highly parallel version due to Bucker & Sauren 1996), 

• transpose-free quasi-minimal residual, TFQMR (Freund 1992), 

• the Chebyshev acceleration, CHEBYSHEV (Young 1981). 

The source code for these methods is distributed with DDSCAT but only PBCGST and PETRKP can 
be called directly via ddscat .par. It is possible (and was done) to add other options by changing 
the code in get f ml . f . A helpful introduction to conjugate gradient methods is provided by the re- 
port "Conjugate Gradient Method Without Agonizing Pain" by Jonathan R. Shewchuk, available as a 
postscript file: |ftp : / /REPORTS .ADM. CS . CMU.EDU/usr0/anon/1994/CMU-CS-94-12 5 .ps| 



In addition to solving the scattering problem for a dipole array, DDSCAT.6.0 can compute the three- 
dimensional force Frad and torque Fiad exerted on this array by the incident and scattered radiation 
fields. This calculation is carried out, after solving the scattering problem, provided DOTORQ has been 
specified in ddscat .par. For each incident polarization mode, the results are given in terms of 
dimensionless efficiency vectors Qpr and Qr, defined by 



15 Calculation of Radiative Force and Torque 



Frad 



(19) 



,pr — 




Qr 



fcTj-ad 



(20) 



■a„ff"rad 
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where Fj-ad and r^ad are the time-averaged force and torque on the dipole array, k = 27r/A is the 
wavenumber in vacuo, and Urad = Eq/8tt is the time-averaged energy density for an incident plane 
wave with amphtude Eq cos{ujt + 0). The radiation pressure efficiency vector can be written 

Qpr = Qcxtk - Qscag , (21) 

where k is the direction of propagation of the incident radiation, and the vector g is the mean direction 
of propagation of the scattered radiation: 

where dH. is the element of solid angle in scattering direction fi, and dCscn/dfl is the differential scat- 
tering cross section. 

Equations for the evaluation of the radiative force and torque are derived by Draine & Weingartner 
(1996). It is important to note that evaluation of Qpr and Qr involves averaging over scattering di- 
rections to evaluate the linear and angular momentum transport by the scattered wave. This evaluation 
requires appropriate choices of the parameters ICTHM and IPHM - see ^22\ 



16 Memory Requirements 

Since Fortran-77 does not allow for dynamic memory allocation, the executable has compiled into it the 
dimensions for a number of arrays; these array dimensions limit the size of the dipole array which the 
code can handle. The source code supplied to you (which can be used to run the sample calculation 
with 192 dipoles) is restricted to problems with targets with a maximum extent of 16 lattice spacings 
along the x-, y-, and z-directions (MXNX=1 6 , MXNY=1 6 , MXNZ=1 6; i.e, the tai-get must fit within 
an 16x 16x 16=4096 cube) and involve at most 9 different dielectric functions (MXC0MP = 9). With this 
dimensioning, the executable requires about 1 .3 MB of memory to run on an Ultrasparc system; memory 
requirements on other hardware and with other compilers should be similar To run larger problems, you 
will need to edit the code to change PARAMETER values and recompile. 

All of the dimensioning takes place in the main program DDSCAT - this should be the only routine 
which it is necessary to recompile. All of the dimensional parameters are set in PARAMETER statements 
appearing before the array declarations. You need simply edit the parameter statements. Remember, of 
course, that the amount of memory allocated by the code when it runs will depend upon these dimen- 
sioning parameters, so do not set them to ridiculously large values! The parameters MXNX, MXNY, MXNZ 
specify the maximum extent of the target in the x-, y-, or z-directions. The parameter MXCOMP speci- 
fies the maximum number of different dielectric functions which the code can handle at one time. The 
comment statements in the code supply all the information you should need to change these parameters. 

The memory requirement for DDSCAT.6.0 (with the netCDF and FFTW options disabled) is ap- 
proximately 

(1690 + 0.764MXNX X MXNY x MXNZ) kbytes (23) 

Thus a 32x32x32 calculation requires 22.5 MBytes, while a 48x48x48 calculation requires 73.0 
MBytes. These values are for the pgf 7 7 compiler on an Intel Xeon system running the Linux op- 
erating system. 



17 Target Orientation 

Recall that we define a "Lab Frame" (LF) in which the incident radiation propagates in the +x direction. 
In dds cat . par one specifies the first polarization state Gqi (which obviously must lie in the y, z plane 
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in the LF); DDSCAT automatically constructs a second polarization state eo2 = x x Gqi orthogonal to 
eoi (here x is the unit vector in the +x direction of the LF. For purposes of discussion we will always let 
unit vectors x, y, z = x x y be the three coordinate axes of the LF. Users will often find it convenient 
to let polarization vectors eoi = y, eo2 — z (although this is not mandatory - see 321> . 

Recall that definition of a target involves specifying two unit vectors, ai and a2, which are imagined 
to be "frozen" into the target. We require a2 to be orthogonal to ai. Therefore we may define a "Target 
Frame" (TF) defined by the three unit vectors ai, 3.2, and 3.3 = ai x a2 . 

For example, when DDSCAT creates a 8x6x4 rectangular solid, it fixes ai to be along the longest 
dimension of the solid, and 3.2 to be along the next-longest dimension. 

Orientation of the target relative to the incident radiation can in principle be determined two ways: 

1 . specifying the direction of ai and a2 in the LF, or 

2. specifying the directions of x (incidence direction) and y in the TF. 

DDSCAT.6.0 uses method 1.: the angles 8, and f3 are specified in the file ddscat . par. The target 
is oriented such that the polar angles 8 and $ specify the direction of ai relative to the incident direction 
x, where the x,y plane has $ = 0. Once the direction of ai is specified, the angle /? than specifies how 
the target is to rotated around the axis ai to fully specify its orientation. A more extended and precise 
explanation follows: 



17.1 Orientation of the Target in the Lab Frame 

DDSCAT uses three angles, 8, <i>, and f3, to specify the directions of unit vectors ai and a2 in the LF 
(see Fig.|6j. 

8 is the angle between ai and x. 

When <i> = 0, ai will lie in the x, y plane. When $ is nonzero, it will refer to the rotation of ai 
around x: e.g., $ = 90° puts ai in the x, z plane. 

When f3 = 0, 3.2 will lie in the x, ai plane, in such a way that when 8 = and $ = 0, a2 is in the 
y direction: e.g, 8 = 90°, $ = 0, /? = has ai = y and a2 ~ — x. Nonzero /? introduces an additional 
rotation of a2 around ai: e.g., 8 = 90°, ^ — 0, f3 = 90° has ai = y and a2 = z. 

Mathematically: 

ai = X cos 8 + y sin 8 cos $ + z sin 8 sin <i> (24) 
£12 = — xsin8cos/3 + y[cos8cos/3cos$ — sin/3sin<i>] 

-|-z [cos 8 cos /? sin <i> + sin (3 cos $] (25) 
as = xsin8sin/3 — y[cos8sin/3cos<l> + cos/3sin<i>] 

— z[cos8sin/3sin<I> — cos/3cos<l>] (26) 



or, equivalently: 



X = ai cos 8 — 32 sin 8 cos f3 + sin 8 sin f3 (27) 
y = ai sin 8 cos $ + a2 [cos 8 cos (3 cos $ — sin /3 sin $] 

— a3[cos8sin/3cos$ + cos/3sin<i>] (28) 
i ~ ki sin 8 sin $ + a2 [cos 8 cos (3 sin $ + sin (3 cos $] 

— 3.3 [cos 8 sin /3 sin $ — cos/? cos $] (29) 



17.2 Orientation of the Incident Beam in the Target Frame 

Under some circumstances, one may wish to specify the target orientation such that x (the direction of 
propagation of the radiation) and y (usually the first polarization direction) and z (= x x y) refer to 
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Figure 6: Target orientation in the Lab Frame, x is the direction of propagation of the incident radiation, and 
y is the direction of the "real" component of the first incident polarization mode. In this coordinate system, 
the orientation of target axis ai is specified by angles and <I>. With target axis ai fixed, the orientation of 
target axis a2 is then determined by angle (3 specifying rotation of the target around ai. When /3 = 0, a2 
hes in the ai,x plane. 
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certain directions in the TF. Given the definitions of the LF and TF above, this is simply an exercise 
in coordinate transformation. For example, one might wish to have the incident radiation propagating 
along the (1,1,1) direction in the TF (example 14 below). Here we provide some selected examples: 

ai, y = a2, z = a3 : © = 0, $ + /3 = 

ai, y = a3, z = -a2 : 6 = 0, $ + /? = 90° 

a2, y = ai, z = -as : e = 90°, /3 = 180°, $ = 

£12, y = aa, z = ai : e = 90°, /3 = 180°, $ = 90° 

as, y = ai, z = 32 : e = 90°, (3 = -90°, $ = 

as, y = a2, z = -ai : e = 90°, (3 = -90°, $ -90° 

-ai, y = aa, z = -as : e = 180°, /3 + $ = 180° 

-ai, y = as, z = a2 : e = 180°, /3 + $ = 90° 

-a2, y = ai, z = as : 6 = 90°, /3 = 0, $ = 

-a2, y = as, z = -ai : e = 90°, /3 = 0, $ = -90° 

-as, y = ai, z = -a2 : © = 90°, (3 = -90°, $ = 

-as, y = a2, z = ai : 6 = 90°, /3 = -90°, $ = 90° 

(ai + a2)/y2, y = a;,, z = (ai - a.2)/^2: 9 = 45°, (3 = 180°, * = 90° 

(ai+a2 + as)/v'3,y= (ai-a2)/V2,z= (ai+a2-2as)/y6: 6 = 54.7356°,/?= 135°, 
= 30°. 

17.3 Sampling in O, and (3 

The present version, DDSCAT.6.0, chooses the angles /3, 9, and $ to sample the intervals (BETAMI , BETAMX), 
(THETMI, THETMX) , (PHIMIN, PHIMAX), where BETAMI, BETAMX, THETMI, THETMX, PHIMIN, 
PHIMAX are specified in ddscat . par . The prescription for choosing the angles is to: 

• uniformly sample in /?; 

• uniformly sample in 

• uniformly sample in cos 9. 

This prescription is appropriate for random orientation of the target, within the specified limits of /3, 
and 9. 

Note that when DDSCAT.6.0 chooses angles it handles (3 and $ differently from 9.^^ The range 
for P is divided into NBETA intervals, and the midpoint of each interval is taken. Thus, if you take 
BETAMI=0, BETAMX=90, NBETA=2 you will get /3 = 22.5° and 67.5°. Similarly, if you take PHIMIN=0, 
PHIMAX=180, NPHI=2 you will get $ = 45° and 135°. 

Sampling in 9 is done quite differently from sampling in j3 and $. First, as already mentioned 
above, DDSCAT.6.0 samples uniformly in cos 9, not 9. Secondly, the sampling depends on whether 
NTHETA is even or odd. 

• If NTHETA is odd, then the values of 9 selected include the extreme values THETMI and THETMX; 
thus, THETMI=0, THETMX=90, NTHETA=3 will give you 9 = 0, 60°, 90°. 

• If NTHETA is even, then the range of cos 9 will be divided into NTHETA intervals, and the mid- 
point of each interval will be taken; thus, THETMI=0, THETMX=90, NTHETA=2 will give you 
9 = 41.41° and 75.52° [cos 9 = 0.25 and 0.75]. 



''This is a change from DDSCAT.4a! !. 
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The reason for this is that if odd NT HE T A is specified, than the "integration" over cos0 is performed 
using Simpson's rule for greater accuracy. If even NTHETA is specified, then the integration over cos Q 
is performed by simply taking the average of the results for the different 9 values. 

If averaging over orientations is desired, it is recommended that the user specify an odd value of 
NTHETA so that Simpson's rule will be employed. 

18 Orientational Averaging 

DDSCAT has been constructed to facilitate the computation of orientational averages. How to go about 
this depends on the distribution of orientations which is apphcable. 

18.1 Randomly-Oriented Targets 

For randomly-oriented targets, we wish to compute the orientational average of a quantity (5(/?, 0, $): 

(Q)^7r^ dp dcosQ d$Q(/3,e,$) . (30) 
Stt^ 7o J-1 Jo 

To compute such averages, all you need to do is edit the file ddscat .par so that DDSCAT knows 
what ranges of the angles /?, 8, and $ are of interest. For a randomly-oriented target with no symmetry, 
you would need to let P run from to 360°, 6 from to 180°, and $ from to 360°. 

For targets with symmetry, on the other hand, the ranges of (3, 8, and $ may be reduced. First of 
all, remember that averaging over $ is relatively "inexpensive", so when in doubt average over to 
360°; most of the computational "cost" is associated with the number of different values of (/3,8) which 
are used. Consider a cube, for example, with axis ai normal to one of the cube faces; for this cube /3 
need run only from to 90°, since the cube has fourfold symmetry for rotations around the axis ai. 
Furthermore, the angle 8 need run only from to 90°, since the orientation (/3,8,$) is indistinguishable 
from (/3, 180° - 8, 360° - $). 

For targets with symmetry, the user is encouraged to test the significance of /3,Q,^ on targets with 
small numbers of dipoles (say, of the order of 100 or so) but having the desired symmetry. 

It is important to remember that DDSCAT .4b handled even and odd values of NTHETA differently 
- see fjsjabove! For averaging over random orientations odd values of NTHETA are to be preferred, as 
this will allow use of Simpson's rule in evaluating the "integral" over cos 8. 

18.2 Nonrandomly- Oriented Targets 

Some special cases (where the target orientation distribution is uniform for rotations around the x axis 
= direction of propagation of the incident radiation), one may be able to use DDSCAT.6.0 with ap- 
propriate choices of input parameters. More generally, however, you will need to modify subroutine 
ORIENT to generate a Ust of NBETA values of /3, NTHETA values of 8, and NPHI values of $, plus two 
weighting arrays WGTA ( 1-NTHETA, 1-NPHI) and WGTB (1-NBETA) . Here WGTA gives the weights 
which should be attached to each (8,$) orientation, and WGTB gives the weight to be attached to each 
/3 orientation. Thus each orientation of the target is to be weighted by the factor WGTAxWGTB. For the 
case of random orientations, DDSCAT.6.0 chooses 8 values which are uniformly spaced in cos 8, and 
P and $ values which are uniformly spaced, and therefore uses uniform weights 
WGTB=1./NBETA 

When NTHETA is even, DDSCAT sets 
WGTA=1./(NTHETAXNPHI) 
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but when NT HE T A is odd, DDSCAT uses Simpson's rule when integrating over 8 and 
WGTA= (1/3 or 4/3 or2/3)/(NTHETAxNPHl) 
Note that the program structure of DDSCAT may not be ideally suited for certain highly oriented 
cases. If, for example, the orientation is such that for a given $ value only one 9 value is possible 
(this situation might describe ice needles oriented with the long axis perpendicular to the vertical in the 
Earth's atmosphere, illuminated by the Sun at other than the zenith) then it is foolish to consider all the 
combinations of Q and $ which the present version of DDSCAT is set up to do. We hope to improve 
this in a future version of DDSCAT. 

19 Target Generation 

DDSCAT contains routines to generate dipole arrays representing targets of various geometries, in- 
cluding spheres, ellipsoids, rectangular solids, cylinders, hexagonal prisms, tetrahedra, two touching 
ellipsoids, and three touching ellipsoids. The target type is specified by variable CSHAPE on line 9 of 
ddscat .par, up to 6 target shape parameters (SHPARl, SHPAR2, SHPAR3, ...) on line 10, and the 
lattice spacing DX D Y D Z on line 1 1 (for a cubic lattice, set DX DY DZ to 1 . 1 . 1 . on hne 11). The 
target geometry is most conveniently described in a coordinate system attached to the target which we 
refer to as the "Target Frame" (TF), so in this section only we will let x,y,z be coordinates in the Target 
Frame. Once the target is generated, the orientation of the target in the Lab Frame is accomplished as 
described in ^171 

Target geometries currently supported include:'^ 

• ELLIPS - Homogeneous, isotropic ellipsoid. 

"Lengths" SHPARl, SHPAR2, SHPAR3 in the x, y, z directions in the TF. (a;/SHPARl)2 + 
(i//SHPAR2)2 + (z/SHPAR3)2 = (f/A, where d is the interdipole spacing. 
The target axes are set to ai = (1, 0, 0) and 3.2 = (0, 1, 0) in the TF. 
User must set NC0MP=1 on line 9 of ddscat .par. 

A homogeneous, isotropic sphere is obtained by setting SHPARl = SHPAR2 = SHPAR3 = 
diameter/d. 

• AN I ELL - Homogeneous, anisotropic elUpsoid. 

SHPARl, SHPAR2, SHPAR3 have same meaning as for ELLIPS. Target axes ai = (1, 0,0) and 
3.2 = (0, 1, 0) in the TF. It is assumed that the dielectric tensor is diagonal in the TF. User must 
set NC0MP=3 and provide xx, yy, zz elements of the dielectric tensor. 

• CONE LL - Two concentric elUpsoids. 

SHPARl, SHPAR2, SHPAR3 specify the lengths along x, y, z axes (in the TF) of the outer 
ellipsoid; SHPAR4, SHPAR5, SHPAR6 ai-e the lengths, along the x, y, z axes (in the TF) of the 
inner ellipsoid. The "core" within the inner ellipsoid is composed of isotropic material 1; the 
"mantle" between inner and outer ellipsoids is composed of isotropic material 2. Target axes 
ai = (1, 0, 0), 3.2 ~ (0, 1, 0) in TF. User must set NC0MP=2 and provide dielectric functions for 
"core" and "mantle" materials. 

• CYLNDR - Homogeneous, isotropic cylinder. 

Length/d=SHPARl, diameter/rf=SHPAR2, with cylinder axis = ai = (1, 0, 0) and a2 = (0, 1,0) 
in the TF. User must set NC0MP=1. 

"in DDSCAT.6.0 not all of the target generation routines have been modified for use with a noncubic lattice. Those routines 
which have not yet been modified (TAR2EL, TAR2SP, TAR3EL, TARCEL, TARHEX, TARTET include code to ensure that they 
are not used with a noncubic lattice; if called with a noncubic lattice, an error message is generated (e.g. "tarhex does not support 
noncubic lattice") and execution is terminated. 
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UNICYL - Homogeneous cylinder with unixial anisotropic dielectric tensor. 

SHPARl, SHPAR2 have same meaning as for CYLNDR. Cylinder axis = ai = (1,0,0), a2 = 
(0, 1, 0). It is assumed that the dielectric tensor e is diagonal in the TF, with Cyy = ezz- User 
must set NC0MP=2. Dielectric function 1 is for E || ai (cylinder axis), dielectric function 2 is for 

E _L ai. 

RCTNGL - Homogeneous, isotropic, rectangular solid. 

X, y, z lengths/d = SHPARl, SHPAR2, SHPAR3. Target axes ai = (1, 0, 0) and a2 = (0, 1, 0) in 
the TF. User must set NC0MP=1. 

HEXGON - Homogeneous, isotropic hexagonal prism. 

Length/d = SHPARl, hexagon side/d = SHPAR2. Target axis ai = (1, 0, 0) in the TF is along the 
prism axis, and target axis a2 = (0,1, 0) in the TF is normal to one of the rectangular faces of the 
hexagonal prism. User must set NC0MP=1. 

TETRAH - Homogeneous, isotropic tetrahedron. 

SHPARl=length/rf of one edge. Orientation: one face parallel to y,z plane (in the TF), opposite 
"vertex" is in +x direction, and one edge is parallel to z axis (in the TF). Target axes ai = (1, 0, 0) 
[emerging from one vertex] and a2 = (0, 1, 0) [emerging from an edge] in the TF. User must set 
NC0MP=1. 

TWOSPH - Two touching homogeneous, isotropic spheroids, with distinct compositions. 

First spheroid has length SHPARl along symmetry axis, diameter SHPAR2 perpendicular to sym- 
metry axis. Second spheroid has length SHPAR3 along symmetry axis, diameter SHPAR4 per- 
pendicular to synmietry axis. Contact point is on Une connecting centroids. Line connecting 
centroids is in x direction. Symmetry axis of first spheroid is in y direction. Symmetry axis of 
second spheroid is in direction y cos(SHPAR5) + z sin(SHPAR5), where y and z are basis vectors 
in TF, and SHPAR5 is in degrees. If SHPAR6=0., then target axes ai = (1, 0, 0), a2 = (0, 1, 0). 
If SHPAR6=1., then axes ai and 3.2 are set to principal axes with largest and 2nd largest moments 
of inertia assuming spheroids to be of uniform density. User must set NC0MP=2 and provide 
dielectric function files for each spheroid. 

TWOELL - Two touching, homogeneous, isotropic ellipsoids, with distinct compositions. 

SHPARl, SHPAR2, SHPAR3=x-length/d, y-length/d, z-\enphld of one elhpsoid. The two el- 
Upsoids have identical shape, size, and orientation, but distinct dielectric functions. The line 
connecting ellipsoid centers is along the x-axis in the TF. Target axes ai = (1, 0, 0) [along line 
cormecting elhpsoids] and a2 = (0, 1, 0). User must set NC0MP=2 and provide dielectric function 
file names for both elhpsoids. Ellipsoids are in order of increasing x: first dielectric function is for 
elhpsoid with center at negative x, second dielectric function for elhpsoid with center at positive 

X. 

TWOAEL - Two touching, homogeneous, anisotropic ellipsoids, with distinct compositions. 

Geometry as for TWOELL; SHPARl, SHPAR2, SHPAR3 have same meanings as for TWOELL. 
Target axes ai = (1, 0, 0) and 3.2 = (0, 1, 0) in the TF. It is assumed that (for both ellipsoids) the 
dielectric tensor is diagonal in the TF. User must set NC0MP=6 and provide xx, yy, zz components 
of dielectric tensor for first ellipsoid, and xx, yy, zz components of dielectric tensor for second 
elhpsoid (ellipsoids are in order of increasing x). 

THRELL - Three touching homogeneous, isotropic ellipsoids of equal size and orientation, 
but distinct compositions. 

SHPARl, SHPAR2, SHPAR3 have same meaning as for TWOELL. Line connecting ellipsoid cen- 
ters is parallel to x axis. Target axis ai = (1, 0, 0) along line of ellipsoid centers, a2 = (0, 1, 0). 
User must set NC0MP=3 and provide (isotropic) dielectric functions for first, second, and third 
elhpsoid. 
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THRAEL - Three touching homogeneous, anisotropic ellipsoids with same size and orienta- 
tion but distinct dielectric tensors. 

SHPARl, SHPAR2, SHPAR3 have same meanings as for THRELL. Target axis ai = (1,0,0) 
along line of ellipsoid centers, a2 = (0, 1, 0). It is assumed that dielectric tensors are all diagonal 
in the TF. User must set NC0MP=9 and provide xx, yy, zz elements of dielectric tensor for first 
ellipsoid, xx, yy, zz elements for second ellipsoid, and xx, yy, zz elements for third elUpsoid 
(ellipsoids are in order of increasing x). 

BLOCKS - Homogeneous target constructed from cubic "blocks". 

Number and location of blocks are specified in separate file blocks, par with following struc- 
ture: 

one line of comments (may be blank) 

PRIN (= or 1 - see below) 

N (= number of blocks) 

B (= width/rf of one block) 

X y z {= position of 1st block in units of Bd) 

xy z{= position of 2nd block) 

X y z {= position of Nth block) 
If PRIN=0, then ai = (1, 0, 0), a2 = (0, 1, 0). If PRIN=1, then ai and a2 are set to principal 
axes with largest and second largest moments of inertia, assuming target to be of uniform density. 
User must set NC0MP=1. 

DWl 9 9 6 - 13 block target used by Draine & Weingartner (1996). 

Single, isotropic material. Target geometry was used in study by Draine & Weingartner (1996) 
of radiative torques on irregular grains, ai and a2 are principal axes with largest and second- 
largest moments of inertia. User must set NC0MP=1 . Target size is controlled by shape parameter 
SHPAR { 1 ) = width of one block in lattice units. 

NSPHER - Multisphere target consisting of the union of N spheres of single isotropic mate- 
rial. 

Spheres may overlap if desired. The relative locations and sizes of these spheres are defined in an 

external file, whose name (enclosed in quotes) is passed through ddscat . par. The length of 

the file name should not exceed 13 characters. The pertinent line in ddscat .par should read 

SHPAR (1) SHPAR (2) '^^Zename ' (quotes must be used) 

where SHPAR { 1 ) = target diameter in x direction (in Target Frame) in units of d 

SHPAR ( 2 ) = to have oi = (1, 0, 0), a2 = (0, 1, 0) in Target Frame. 

SHPAR ( 2 ) = 1 to use principal axes of moment of inertia tensor for ai (largest I) and 02 (inter- 
mediate /). 

filename is the name of the file specifying the locations and relative sizes of the spheres. 
The overall size of the multisphere target (in terms of numbers of dipoles) is determined by pa- 
rameter SHPAR ( 1 ) , which is the extent of the multisphere target in the x-direction, in units of 
the lattice spacing d. The file 'filename ' should have the following structure: 

N (= number of spheres) 

one fine of comments (may be blank) 

xi yi zi ai (arb. units) 

X2 t/2 -22 (i2 (arb. units) 

xn Vn zn (In (arb. units) 
where Xj, yj, Zj are the coordinates of the center, and aj is the radius of sphere j. 
Note that Xj, yj, Zj, aj (J = 1,...,N) estabhsh only the shape of the iV— sphere target. For in- 
stance, a target consisting of two touching spheres with the Une between centers parallel to the x 
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axis could equally well be described by lines 3 and 4 being 

0.5 

1 0.5 

or 

000 1 
200 1 

The actual size (in physical units) is set by the value of Ues specified in ddscat . par, where, as 
always, oieff = (3y^/47r)^/^, where V is the volume of material in the target. 
User must set NC0MP=1. 

PRISM3 - Triangular prism of homogeneous, isotropic material. 

SHPARl, SHPAR2, SHPAR3, SHPAR4 = a/rf, &/a, c/a, L/a 

The triangular cross section has sides of width a, b, c. L is the length of the prism, d is the lattice 
spacing. The triangular cross-section has interior angles a, 0, 7 (opposite sides a, b, c) given by 
cos a = (&2 + - a^)/2bc, cos ,'3 = (a^ + c? - b'^)/2ac, cos 7 = (a^ + 6^ _ c?)/2ab. In the 
Target Frame, the prism axis is in the x direction, the normal to the rectangular face of width a is 
(0,1,0), the normal to the rectangular face of width b is (0, — C0S7, sin 7), and the normal to the 
rectangular face of width c is (0, — cos 7, — sin 7). 

LYRSLB - Multilayer rectangular slab with overall x, y, z lengths = SHPARl xd, ay = SHPAR2 x d, 
Uz = SHPARSxd, with the boundary between layers 1 and 2 occuring at xi — SHPAR4 x Ux, 
boundary between layers 2 and 3 occuring at X2 = SHPAR5 x Ux [if SHPAR5> 0], and boundary 
between layers 3 and 4 occuring at xz = SHPAR6 x Ux [if SHPAR6> 0], To create a simple bi- 
layer slab, just set SHPAR5=SHPAR6=0. User must set NCOMP=2,3, or 4 and provide dielectric 
function files for each of the two layers. 

FRMFIL - List of dipole locations and "compositions" obtained from a file. 

This option causes DDSCAT to read the target geometry information from a file shape . dat 
instead of automatically generating one of the geometries listed above. The shape . dat file 
is read by routine REASHP (file reashp. f). The user can customize REASHP as needed to 
conform to the manner in which the target geometry is stored in file shape . dat. However, as 
supplied, REASHP expects the file shape . dat to have the following structure: 

- one Une containing a description; the first 67 characters wiU be read and printed in various 
output statements. 

- N = number of dipoles in target 

- ciix ci\y a\z = x,y,z components of ai 

- 0'2x 0'2y 0'2z = x,y,z Components of a2 

- (line containing comments) 

- dummy IXYZ (1, 1) IXYZ(1,2) IXYZ(1,3) IC0MP(1,1) IC0MP(1,2) IC0MP(1,3) 

- dummy IXYZ {2,1) IXYZ(2,2) IXYZ(2,3) IC0MP(2,1) ICOMP(2,2) ICOMP(2,3) 

- dummy IXYZ (3,1) IXYZ (3, 2) IXYZ (3, 3) IC0MP(3,1) ICOMP(3,2) ICOMP(3,3) 

- dummy IXYZ {J, 1) IXYZ (J, 2) IXYZ (J, 3) ICOMP(J, 1) ICOMP(J, 2) ICOMP(J, 3) 



- dummy IXYZ (N, 1) IXYZ(N,2) IXYZ(N,3) IC0MP(N,1) IC0MP(N,2) IC0MP(N,3) 

The user should be able to easily modify these routines, or write new routines, to generate targets 
with other geometries. The user should first examine the routine target . f and modify it to call any 
new target generation routines desired. Alternatively, targets may be generated separately, and the target 
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description (locations of dipoles and "composition" corresponding to x,y,z dielectric properties at each 
dipole site) read in from a file by invoking the option FRMFIL in ddscat . f . 

It is often desirable to be able to run the target generation routines without running the entire 
DDSCAT code. We have therefore provided a program CALL TARGET which allows the user to gener- 
ate targets interactively; to create this executable just type'^ 
make calltarget . 

The program calltarget is to be run interactively; the prompts are self-explanatory. You may need 
to edit the code to change the device number IDVOUT as for DDSCAT (see i|6.1l above). 

After running, calltarget will leave behind an ASCII file target . out which is a list of the 
occupied lattice sites in the last target generated. The format of target . out is the same as the format 
of the shape . dat files read if option FRMFIL is used (see above). Therefore you can simply 

mv target. out shape.dat 
and then use DDSCAT with the option FRMFIL in order to input a target shape generated by CALLTARGET. 



20 Scattering Directions 

DDSCAT calculates scattering in selected directions, and elements of the scattering matrix are reported 
in the output files vxxryykzzz . sea . The scattering direction is specified through angles 6s and 0^ (not 
to be confused with the angles 8 and $ which specify the orientation of the target relative to the incident 
radiation!). 

The scattering angle Bg is simply the angle between the incident beam (along direction x) and the 
scattered beam (dg — for forward scattering, 9s — 180° forbackscattering). 

The scattering angle (ps specifies the orientation of the "scattering plane" relative to the ic, y plane 
in the Lab Frame. When (j>s = the scattering plane is assumed to coincide with the x, y plane. When 
cf>s = 90° the scattering plane is assumed to coincide with the x, z plane. Within the scattering plane the 
scattering directions are specified by < 6*^ < 180°. 

Scattering directions for which the scattering properties are to be calculated are set in the file 
ddscat .par by specifying one or more scattering planes (determined by the value of (ps) and for 
each scattering plane, the number and range of 6s values. The only limitation is that the number of 
scattering directions not exceed the parameter MXSCA in DDSCAT . f (in the code as distributed it is set 
toMXSCA=1000). 



21 Incident Polarization State 

Recall that the "Lab Frame" is defined such that the incident radiation is propagating along the x axis. 
DDSCAT.6.0 allows the user to specify a general elliptical polarization for the incident radiation, by 
specifying the (complex) polarization vector Gqi. The orthonormal polarization state eo2 = x x Gqi is 
generated automatically if ddscat . par specifies I0RTH=2. 

For incident linear polarization, one can simply set Gqi = y by specifying (0,0) (1,0) (0,0) 
in dds cat .par; then eo2 = z For polarization mode Gqi to correspond to right-handed circular polar- 
ization, set eoi = (y+iz)/V2 by specifying (0,0) (1,0) (0,1) in ddscat .par (DDSCAT.6.0 
automatically takes care of the normalization of eoi); then eo2 — {~iy + z) /y/2, corresponding to left- 
handed circular polarization. 

"Non-Unix sites: The source code for CALLTARGET is in the file MISC . FOR. You must compile and link this to ERRMSG, 
GETSET, REASHP, TAR2EL, TAR2SP, TAR3EL, TARBLOCKS, TARCEL, TARCYL, TARELL, TARGET, TARHEX, TARREC, 
TARTET, andWRIMSG. The source code for ERRMSG is in SRC2 . FOR, GETSET is in SRC5 . FOR, and the rest are in SRC8 . FOR 
and SRC9 . FOR . 
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22 Averaging over Scattering: g{l) = (cos Og), etc. 



DDSCAT automatically carries out numerical integration of various scattering properties. In particular, 
it calculates the mean value g(l) = (cos^s) for the scattered intensity for each incident polarization 
state. This is accomplished by evaluating the scattered intensity for ICTHM different values of 9s (in- 
cluding 9s — and 9s = n), and taking a weighted sum. For each value of 9s except and tt, the 
scattering intensity is evaluated for IPHM different values of the scattering angle (f>s. The angular in- 
tegration over 9s is accomplished using Simpson's rule (with uniform intervals in cos9s), so ICTHM 
should be an odd number The angular integration over (j)s is accomplished by sampling uniformly in (j)s 
with uniform weighting; IPHM can be either even or odd. 

The following quantities are evaluated by this angular integration: 

• g = (cos 6's)x + (sin 9s cos (f>s)y + (sin 9s sin (f>s)z (see ^151 : 

• Qr(seefIT5}. 

It is important that the user recognize that accurate evaluation of these angular averages requires ade- 
quate sampling over scattering angles. For small values of the size parameter x = 27raefF / A, the angular 
distribution of scattered radiation has a dipolar character and the sampling in 9s and (j)s does not need to 
be very fine, so ICTHM and IPHM need not be large. For larger values of the size parameter x, however, 
higher multipoles in the scattered radiation field become important, and finer sampling in 9s and (j>s is 
required. We do not have any foolproof prescription to offer, since the scattering pattern will depend 
upon the target geometry and dielectric constant in addition to overall size parameter. However, as a 
very rough guide, we suggest that the user specify values of ICTHM and IPHM satisfying 



The sample ddscat .par file supplied has ICTHM = 33 and IPHM = 12; the above criteria would 
suggest that this would be suitable for x < 5. 

The cpu time required for evaluation of these angular averages is proportional to [2 + IPHM(ICTHM — 
2)]. Since the computational time spent in evaluating these angular integrals can be a significant part 
of the total, it is important to choose values of ICTHM and IPHM which will provide a suitable balance 
between accuracy (in this part of the overall calculation) and cpu time. 

Within one scattering plane, the scattered intensity tends to have approximately {1 + x) peaks for 
< 9s < TT, so that the above prescription for ICTHM would have at least 5 sampling points per 
maximum. The angular distribution over (ps is usually not as structured as that over 9s so we suggest 
that IPHM need not be as large as ICTHM. We have refrained from "hard-wiring" the values of ICTHM 
and IPHM because we are not confident of the reliability of the recommended criteria ( I31I32> - it is 
up to the user to specify appropriate values of ICTHM and IPHM according to the requirements of the 
problem being addressed. 

23 Mueller Matrix for Scattering in Selected Directions 
23.1 Two Orthogonal Incident Polarizations (I0RTH=2) 

DDSCAT.6.0 internally computes the scattering properties of the dipole array in terms of a complex 
scattering matrix fmi{Ss, <t>s) (Draine 1988), where index I = 1,2 denotes the incident polarization state, 
m = 1, 2 denotes the scattered polarization state, and 9s, (ps specify the scattering direction. Normally 
DDSCAT is used with I0RTH=2 in ddscat .par, so that the scattering problem will be solved for 
both incident polarization states (I — 1 and 2); in this subsection it will be assumed that this is the case. 



ICTHM > 5(1 + a;) 
IPHM > 2(1 + a;) 



(31) 
(32) 
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Incident polarization states 1 — 1,2 correspond to polarization states eoi, 692; recall that polarization 
state eoi is user-specified, and eo2 = x x Bq^. Scattered polarization state to = 1 corresponds to 
linear polarization of the scattered wave parallel to the scattering plane (ei — ep — 9s) and m — 2 
corresponds to linear polarization perpendicular to the scattering plane (in the +0s direction). The 
scattering matrix was defined (Draine 1988) so that the scattered electric field is related to the 
incident electric field (0) at the origin (where the target is assumed to be located) by 



E^ • \ _ exp(ik • r) / /ii /12 \ / E,(0) • eoi 



Es • 0s / kr V /21 /22 / V • ^02 



(33) 



The 2x2 complex amplitude scattering matrix (with elements Si, S2, S3, and S4) is defined so that (see 
Bohren & Huffman 1983) 



E, • \ exp(ik ■ r) r S2 S3\ f E,(0) • e, 
~Es-4>s J -ikr \ Si Si )\ E,(0)-e,_ 



(34) 



where j| , are (real) unit vectors for incident polarization parallel and perpendicular to the scattering 
plane (with the customary definition of e^^ — e.^\ x x). 



From ( I33I34> we may write 



S2 ^3 W E,(0)-e,|| \ _ / W E,(0) 

^4 ^1 M E,(0)-e,i i M -/21 -/22 / I E,(0) 



=01 
-02 



(35) 



Let 



a = 


^01 


■y , 


(36) 


h = 


^01 




(37) 


c = 


^02 


■y , 


(38) 


d = 


p* 
'='02 


• z . 


(39) 



Note that since eoi, eo2 could be complex (i.e., elUptical polarization), the quantities a, b, c, d are com 
plex. Then 

c d I \ z 



(40) 



and eq. J35> can be written 

^2 ^3 W E,(0) • e,|| \ / -fii -/12 \f a 6 \ / E,(0) • y 
^4 Si ){E,{0)-e,^ J \ hi /22 )\c d){-E,{0)-i 

The incident polarization states ei|| and e^j^ are related to y, z by 

y \ ^ ( cos 0s sin 0s \ f e^ii 
z J \ sin 0s -cos 0s J \ eij_ 



(41) 



(42) 



substituting J42t into ( I41> we obtain 

f S2 S3\f E,(0) • e,t| \ . / -fii -fi2 \f a 6 \ / cos 0s sin 0s \ ( E,(0) • e,,, ^ 
\Si Si J \ E,(0) -e,^ J \ f2i f22 )\c d )\ sin0s -cos0s ) \ E,(0) • e,^ ^ 

(43) 
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Eq. ( I43> must be true for all Ei(0), so we obtain an expression for the complex scattering amplitude 
matrix in terms of the fmi : 

( ? ?')=^fV" "'"t ) ■ (44) 

\ S4 Si J V /21 h2 ) \ c d J \ sm0s -cos 0s J 

This provides the 4 equations used in subroutine GETMUELLER to compute the amplitude scattering 
matrix elements: 

51 = [/2i(6cos0s - asin^s) + /22(dcos0s - csin^s)] , (45) 

52 = [/ii(acos(?!)s + fosin^s) + /i2(ccos(?!)5 + dsin^s)] , (46) 
S-i = i [/ii(6cos0s - asin^s) + /i2(rfcos0s - csin^s)] , (47) 
Sa = i [/2i(acos0s + fosin^s) + /22(ccos0s + dsin^s)] . (48) 



23.2 Stokes Parameters 

It is both convenient and customary to characterize both incident and scattered radiation by 4 "Stokes 
parameters" - the elements of the "Stokes vector". There are different conventions in the literature; we 
adhere to the definitions of the Stokes vector (I,Q,U,V) adopted in the excellent treatise by Bohren & 
Huffman (1983), to which the reader is referred for further detail. Here are some examples of Stokes 
vectors (/, Q, [/, V) = (1, Q/I, U/I, V/I)I: 

• (1,0, 0,0)/: unpolarized light (with intensity /); 

• (1,1,0, 0)/ : 100% linearly polarized with E parallel to the scattering plane; 

• (1, —1, 0, 0)/ : 100% linearly polarized with E perpendicular to the scattering plane; 

• (1, 0, 1, 0)/ : 100% linearly polarized with E at +45° relative to the scattering plane; 

• (1, 0, —1, 0)/ : 100% linearly polarized with E at -45° relative to the scattering plane; 

• (1,0, 0, 1)/ : 100% right circular polarization {i.e., negative helicity); 

• (1, 0, 0, —1)/ : 100% left circular polarization {i.e., positive helicity). 



23.3 Relation Between Stokes Parameters of Incident and Scattered Radiation: 
The Mueller Matrix 

It is convenient to describe the scattering properties in terms of the 4x4 Mueller matrix relating the 
Stokes parameters {Ii,Qi,Ui,Vi) and {Ig, Qs,Us,Vs) of the incident and scattered radiation: 





( 


Qs 


1 


Us 






V 



Su S12 5*13 5*14 

521 5*22 '5'23 5*24 

S31 5*32 5*33 5*34 

y 541 5*42 5*43 5*44 



\ 


( \ 




Q^ 




u, 







(49) 



Once the amplitude scattering matrix elements are obtained, the Mueller matrix elements can be com- 
puted (Bohren & Huffman 1983): 



Sii 
S12 

Sl3 
Sl4 



= {\Si 



|52p + |53p 



54'"^ 



(|S'2p - IS*!]' -I- |04 
Re (52^3* +5154*) 

Im(5253*-5i54*) 



l^4p)/2 , 
l^3p)/2 , 
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<5'21 


- (.\S2f-\S^ 


\' + \S3f 


S22 


= i\Si\' + \S, 


\'-\S3\' 


823, 


= Re {S2S; ~ 


Si SI) , 


824, 


= Im(52^3 + 


^1^4*) , 


S31 


= Re(S'2S'4* + 


SiS*,) , 


S32 


= Re(52^4*- 


Sis;) , 


S33 


= Re(S'iS'* + 


S3S4) , 


S34 


= Im {S2SI + 


S4S2) , 


S41 


^ Im {S4S; + 


^1^3) , 


S42 


— Im (<S'4S'2 — 


^1^3) , 


S43 


= Im (S'iS'2 — 


^3^4*) , 


S44 


= Re (S'iS'2* - 


S3S4) . 



|S4|^, 
|S4p- 



/2 
/2 



(50) 

These matrix elements are computed in DDSCAT and passed to subroutine WRITESCA which handles 
output of scattering properties. Although the Muller matrix has 16 elements, only 9 are independent. 

The user can select up to 9 distinct Muller matrix elements to be printed out in the output files 
wxxryykzzz . sea and vixxryyori . avg; this choice is made by providing a list of indices in ddscat . par 
(see AppendixIXIi. 

If the user does not provide a list of elements, WRITESCA will provide a "default" set of 6 selected 
elements: S'n, S'21, S31, S'41 (these 4 elements describe the intensity and polarization state for scattering 
of unpolarized incident radiation), S'12, and S'13. 

In addition, WRITESCA writes out the linear polarization P of the scattered light for incident unpo- 
larized light: 



P 



Sii 



(51) 



23.4 One Incident Polarization State Only (I ORTH= 1) 

In some cases it may be desirable to limit the calculations to a single incident polarization state - for 
example, when each solution is very time-consuming, and the target is known to have some symmetry 
so that solving for a single incident polarization state may be sufficient for the required purpose. In this 
case, set I0RTH=1 in ddscat .par. 

When I0RTH=1, only /n and /21 are available; hence, DDSCAT cannot automatically generate 
the Mueller matrix elements. In this case, the output routine WRITESCA writes out the quantities 
1/21 p, Re(/ii/2i), and Im(/ii/|4) for each of the scattering directions. 

Note, however, that if IPHI is greater than 1, DDSCAT will automatically set I0RTH=2 even if 
ddscat . par specified I0RTH=1: this is because when more than one value of the target orientation 
angle <I> is required, there is no additional "cost" to solve the scattering problem for the second incident 
polarization state, since when solutions are available for two orthogonal states for some particular target 
orientation, the solution may be obtained for another target orientation differing only in the value of 
$ by appropriate linear combinations of these solutions. Hence we may as well solve the "complete" 
scattering problem so that we can compute the complete Mueller matrix. 
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24 Using MPI (Message Passing Interface) 



Many (probably most) users of DDSCAT will wish to calculate scattering and absorption by a given tar- 
get for more than one orientation relative to the incident radiation (e.g., when calculating orientational 
averages). DDSCAT can automatically carry out the calculations over different target orientations. Nor- 
mally, these calculations will be carried out sequentially. However, on a multiprocessor system [either 
a symmetric multiprocessor (SMP) system, or a cluster of networked computers] with MP I (Message 
Passing Interface) software installed, DDSCAT.6.0 allows the user to carry out parallel calculations of 
scattering by different target orientations, with the results gathered together and with orientational aver- 
ages calculated just as they are when sequential calculations are carried out. Note that the MPI-capable 
executable can also be used for ordinary serial calculations using a single cpu. 

More than one implementation of MPI exists. MPI support within DDSCAT.6.0 is compliant with 
MPI- 1 .2 and MPI-2 standards and should be usable under any implementation of MP I that is compat- 
ible with those standards. 

At Princeton University Department of Astrophysical Sciences we are using MPICh'^, a publicly- 
available implementation of MP I. 



24.1 Source Code for the MPI-compatible version of DDSCAT 

In order to use DDSCAT with MPI, the user must first make sure that the correct version of the code is 
being used. To prepare the MPI-capable version: 

1. In DDSCAT . f : 

• Make sure the line 

INCLUDE 'mpif.h' 
is not commented out. 

• Make sure the line 

C INTEGER MPI_COMM_WORLD 

is commented out. 

2. Inmpisubs.f, SUBROUTINE COLSUM : 

• Make sure the line 

INCLUDE 'mpif.h' 
is not commented out. 

• Make sure the line 

C INTEGER MPI_COMM_WORLD,MPI_SUM,MPI_REAL,MPI_COMPLEX 

is commented out. 

3. Inmpisubs.f, SUBROUTINE SHAREl : 

• Make sure the line 

INCLUDE 'mpif.h' 
is not commented out. 

• Make sure the lines 

C INTEGER MPI_CHARACTER,MPI_COMM_WORLD,MPI_INTEGER, 

C & MPI_INTEGER2,MPI_REAL,MPI_C0MPLEX 

are commented out. 

4. Inmpisubs.f, SUBROUTINE SHARE2 : 

http : / /www . mpi-f orum. org/| 

http : / / www . mcs . ani . gov/mpi/mpich| 
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• Make sure the line 

INCLUDE 'mpif.h' 
is not commented out. 

• Make sure the line 

C INTEGER MPI_COMM_WORLD,MPI_REAL,MPI_COMPLEX 

is commented out. 

24.2 Compiling and Linking the MPI-compatible version of DDSCAT 

In the Makefile, enable the following definitions: 

MPIsrc = mpi.subs . f 

MPIobj = mpi_subs.o 

At Princeton University Dept. of Astrophysical Sciences we are currently using mpich-1 . 2 . 4^". 
mpich provides a "compiler wrapper" mpif 7 7 that is supposed to automatically set the variables 
MP I CH_F 7 7 and MP I CH_F 7 7 L INKER to values appropriate for use of the pg f 7 7 compiler from PG? ' 
which is installed on our beowulf cluster Other compilers, e.g., the Intel f77 compiler, can also be em- 
ployed. 

24.3 Code Execution Under MPI 

Local installations of MPI will vary - you should consult with someone familiar with way MPI is 
installed and used on your system. 

At Princeton University Dept. of Astrophysical Sciences we use PBS (Portable Batch System)^^ to 
schedule jobs. MPI jobs are submitted using PBS by first creating a shell script such as the following 
example file pbs . submit: 

# ! /bin/bash 

#PBS -1 nodes=2 :ppn=l 

#PBS -1 mem=12 00MB, pmem=300MB 

#PBS -m bea 

#PBS -j oe 

Cd $PBS_0_WORKDIR 

/usr/local/bin/mpiexec ddscat 

The lines beginning with #PBS -1 specify the required resources: 

#PBS -1 nodes=2 : ppn=l specifies that 2 nodes are to be used, with 1 processor per node. 
#PBS -1 mem=l 2 0MB, pmem=3 0MB specifies that the total memory required (mem) is 1200MB, 
and the maximum physical memory used by any single process (pmem) is 300MB. The actual definition 
of mem is not clear, but in practice it seems that it should be set equal to 2 x (node s) x (ppn) x (pmem) . 
#PBS -m bea specifies that PBS should send email when the job begins (b), and when it ends (e) or 
aborts (a). 

#PBS - j oe specifies that the output from stdout and stderr will be merged, intermixed, as stdout. 
This example assumes that the executable dscat is located in the same directory where the code is to 
execute and write its output. If ddscat is located in another directory, simply give the full pathname, 
to it. The qsub command is used to submit the PBS job: 

% qsub pbs. submit 

I http : / /www-unix .mcs . ani . gov/mpi/mpich/| 
http : / /www. pgroup . com| 
^ jhttp : / /www . openpbs . org| 
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As the calculation proceeds, the usual output files will be written to this directory: for each wave- 
length, target size, and target orientation, there will be a file waarbbkccc . sea, where aa=00, 01, 02, 
... specifies the wavelength, bb=00, 01, 02, ... specifies the target size, and cc=000, 001, 002, ... spec- 
ifies the orientation. For each wavelength and target size there will also be a file waarbbori . avg 
with orientationally-averaged quantities. Finally, there will also be tables qtable, and qtable2 with 
orientationally-averaged cross sections for each wavelength and target size. 

In addition, each processor employed will write to its own log file ddscat . logjinn, where 
nnn=000, 001, 002, .... These files contain information concerning cpu time consumed by different 
parts of the calculation, convergence to the specified error tolerance, etc. If you are uncertain about how 
the calculation proceeded, examination of these log files is recommended. 

25 Graphics and Postprocessing 

25.1 IDL 

At present, we do not offer a comprehensive package for DDSCAT data postprocessing and graphical 
display in IDL. However, there are several developments worth mentioning: First, we offer several out- 
put capabilities from within DDSCAT: ASCII (see ^ToTl . FORTRAN unformatted binaiy (see JTa2b . 
and netCDF portable binary (see i|10.3t . Second, we offer several skeleton IDL utilities: 

• bhmie .pro is our translation to IDL of the popular Bohren-Huffman code which calculates 
efficiencies for spherical particles using Mie theory. 

• readbin.pro reads the FORTRAN unformatted binary file written by routine writebin.f. 
The variables are stored in a common block. 

• readnet.pro reads NetCDF portable binary file and should be the method of choice for IDL 
users. It offers random data access. 

• mie . pro is an example of an interface to binary files and the Bohren-Huffman code. It plots a 
comparison of DDSCAT results with scattering by equivalent radius spheres. 

At present the IDL code is experimental. 

26 Miscellanea 

Additional source code, refractive index files, etc., contributed by users will be located in the directory 
DDA/misc. These routines and files should be considered to be not supported by Draine and Flatau - 
caveat receptor! These routines and files should be accompanied by enough information (e.g., comments 
in source code) to explain their use. 

27 Finale 

This User Guide is somewhat inelegant, but we hope that it will prove useful. The structure of the 
ddscat .par file is intended to be simple and suggestive so that, after reading the above notes once, 
the user may not have to refer to them again. 

The file rel_notes in DDA/ doc lists known bugs in DDSCAT.6.0 . Up-to-date release notes will 
be available at the DDSCAT web site, 

|http : //www .astro . princeton . edu/ ^draine /DP SCAT . html 

Users are encouraged to provide B. T. Draine (draine@astro . princeton . edu) with their 
email address; bug reports, and any new releases of DDSCAT, will be made known to those who do! 
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p. J. Flatau maintains the WWW page "SCATTERLIB - Light Scattering Codes Library" with URL 
fhttp ://atol.ucsd . edu/~pf latau The SCATTERLIB Internet site is a library of light scat- 
tering codes. Emphasis is on providing source codes (mostly FORTRAN). However, other information 
related to scattering on spherical and non-spherical particles is collected: an extensive list of references 
to light scattering methods, refractive index, etc. This URL page contains section on the discrete dipole 
approximation. 

Concrete suggestions for improving DDSCAT (and this User Guide) are welcomed. If you wish to 
cite this User Guide, we suggest the following citation: 

Draine, B.T., & Flatau, RJ. 2000, "User Guide for the Discrete Dipole Approximation Code DDSCAT 
(Version 6.0)", http : //a rxi v. org/ ab s/ as tro-ph/? ? ? ? ? ? 

Finally, the authors have one special request: We would very much appreciate preprints and (espe- 
cially!) reprints of any papers which make use of DDSCAT! 
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A Understanding and Modifying dd scat. par 



In order to use DDSCAT to perform the specific calculations of interest to you, it will be necessary to 
modify the ddscat . par file. Here we list the sample ddscat . par file, followed by a discussion of 
how to modify this file as needed. Note that all numerical input data in DDSCAT is read with free-format 
READ{IDEV, *) . . . statements. Therefore you do not need to worry about the precise format in which 
integer or floating point numbers are entered on a line. The crucial thing is that lines in ddscat . par 
containing numerical data have the correct number of data entries, with any informational coimnents 
appearing after the numerical data on a given Une. 

' =================== Parameter file ===================' 

'**** PRELIMINARIES ****' 

'NOTORQ'= CMT0RQ*6 (DOTORQ, NOTORQ) — either do or skip torque calculations 
'PBCGST'= CMDS0L*6 (PBCGST, PETRKP) — select solution method 
'GPFAFT'= CMETHD*6 (GPFAFT, FFTW21, CONVEX) 

'LATTDR'= CALPHA*6 (LATTDR is only supported option now) 
'NOTBIN'= CBINFLAG (ALLBIN, ORIBIN, NOTBIN) 
'NOTCDF'= CNETFLAG (ALLCDF, ORICDF, NOTCDF) 

' RCTNGL' = CSHAPE*6 (FRMFIL, ELLIPS, CYLNDR, RCTNGL, HEXGON, TETRAH, UNICYL, UNIELL) 
8. 6. 4. = shape parameters PARI, PAR2, PAR3 

1 = NCOMP = number of dielectric materials 

' TABLES' = CDIEL*6 (TABLES, H20ICE, H20LIQ; if TABLES, then filenames follow...) 
'diel.tab' = name of file containing dielectric function 
'**** CONJUGATE GRADIENT DEFINITIONS ****' 

= INIT (TO BEGIN WITH | X0> = 0) 

l.OOe-5 = ERR = MAX ALLOWED (NORM OF | G>=AC | E>-ACA | X> ) / (NORM OF AC|E>) 
'**** ANGLES FOR CALCULATION OF CSCA, G' 

33 = ICTHM (number of theta values for evaluation of Csca and g) 

12 = IPHM (number of phi values for evaluation of Csca and g) 

' **** Wavelengths (micron) 

6.283185 6.283185 1 ' INV = wavelengths (first, last, how many, how=LIN, INV, LOG) 
'**** Effective Radii (micron) **** ' 

1. 1. 1 'LIN' = eff. radii (first, last, how many, how=LIN, INV, LOG) 
' **** Define Incident Polarizations ****' 

(0,0) (l.,0.) (0.,0.) = Polarization state eOl (k along x axis) 

2 = lORTH (=1 to do only pol. state eOl; =2 to also do orth. pol . state) 

1 = IWRKSC (=0 to suppress, =1 to write ".sea" file for each target orient. 
'**** Prescribe Target Rotations ****' 

0. 0. 1 = BETAMI, BETAMX, NBETA (beta=rotation around al) 
0. 90. 3 = THETMI, THETMX, NTHETA (theta=angle between al and k) 
0. 0. 1 = PHIMIN, PHIMAX, NPHI (phi=rotation angle of al around k) 
'**** Specify first IWAV, TRAD, lORI (normally 0) ****' 
= first IWAV, first TRAD, first TORI (0 to begin fresh) 
'**** Select Elements of S_ij Matrix to Print ****' 

6 = NSMELTS = number of elements of S_ij to print (not more than 9) 
11 12 21 22 31 41 = indices ij of elements to print 
' **** Specify Scattered Directions ****' 

0. 0. 180. 30 = phi, thetan_min, thetan_max, dtheta (in degrees) for plane A 
90. 0. 180. 30 = phi, ... for plane B 
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Lines Comments 

1-2 comment lines - no need to change. 

3 NOTORQ if torque calculation is not required; 
DOTORQ if torque calculation is required. 

4 PBCGST is recommended; other option is PETRKP (see fl4l . 

5 GPFAFT is supplied as default, but FFTW2 1 is recommended if DDSCAT has been compiled with 
FFTW support (see ii ii6.5l[T3t : only other option is CONVEX (valid only on a Convex computer). 

6 LATTDR (Draine & Goodman [2000] Lattice Dispersion Relation is only option now supported) 
(see SHi. 

7 ALLBIN for full unformatted binary dump ( i|10.2l l; 

ORIBIN for unformatted binary dump of orientational averages only; 
NOTBIN for no unformatted binary output. 

8 ALLCDF for output in netCDF format (must have netCDF option enabled; cf. ill0.3> : 
ORICDF for orientational averages in netCDF format (must have netCDF option enabled). 
NOTCDF for no output in netCDF format; 

9 specify choice of target shape (see ijlllfor description of options RCTNGL, ELLIPS, TETRAH, ...) 

10 shape parameters SHPARl, SHPAR2, SHPAR3, ... (see 

1 1 number of different dielectric constant tables ( ^121 . 

12 TABLES - dielectric function to be provided via input table. 

13 name(s) of dielectric constant table(s) (one per line). 

14 comment line - no need to change. 

15 is recommended value of parameter INIT. 

16 ERR = eiTor tolerance h: maximum allowed value of \ A^E — AP\/\A^ E\ [see eq.lfTstl. 

17 comment line - no need to change. 

18 ICTHM - number of 6s values for angular averages ( i|22> . 

19 IPHM - number of <j)s values for angular averages. 

20 comment line - no need to change. 

21 A - first, last, how many, how chosen. 

22 comment line - no need to change. 

23 OefF - first, last, how many, how chosen. 

24 comment line - no need to change. 

25 specify x,y,z components of (complex) incident polarization Gqi ( ^2li 

26 lORTH = 1 to do one polarization state only; 

2 to do second (orthogonal) incident polarization as well. 

27 IWRKSC = to suppress writing of ".sea" files; 
2 to enable writing of ".sea" files. 

28 comment line - no need to change. 

29 j3 (see i|17l i - first, last, how many . 

30 Q - first, last, how many. 

31 <& - first, last, how many. 

32 comment line - no need to change. 

33 IWAVO IRADO lORlO - specify Starting values of integers I WAV IRAD I ORI (normally 0). 

34 comment line - no need to change. 

35 Ns = number of scattering matrix elements (must be < 9) 

36 indices ij of Ns elements of the scattering matrix Sij 

37 comment line - no need to change. 

38 (j)s for first scattering plane, 9s^min, Os,max, how many Og values; 
39,... (j)s for 2nd,... scattering plane, ... 
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B wx\;r3ryori . avg Files 



The filewOOrOOori. avg contains the results for the first wavelength (wOO) and first target radius 
(rOO) averaged over orientations (ori . avg). The wOOrOOori. avg file generated by the sample 
calculation should look like the following: 

DDSCAT DDSCAT.6.0 [03.07.13] 

TARGET Rectangular prism; NX,NY,NZ= 8 6 4 

GPFAFT method of solution 

LATTDR prescription for polarizabilies 

RCTNGL shape 

192 = NATO = number of dipoles 



AEFF= 1.00000 = effective radius (physical units) 

WAVE= 6.28319 = wavelength (physical units) 
K*AEFF= 1.00000 = 2*pi*aeff /lambda 

n= ( 1.3300 , 0.0100), eps.= ( 1.7688 , 0.0266) |m|kd= 0.3716 for subs. 1 
TOL= l.OOOE-05 = error tolerance for CCG method 
ICTHM= 33 = theta values used in comp. of Qsca,g 
IPHIM= 12 = phi values used in comp. of Qsca,g 
( 1.00000 0.00000 0.00000) = target axis Al in Target Frame 
( 0.00000 1.00000 0.00000) = target axis A2 in Target Frame 
( 0.27942 0.00000 0.00000) = k vector (latt. units) in Lab Frame 
( 0.00000, 0.00000) ( 1.00000, 0.00000) ( 0.00000, O.00000)=inc.pol.vec. 1 inLF 
( 0.00000, 0.00000) ( 0.00000, 0.00000) ( 1.00000, O.00000)=inc.pol.vec. 2 in LF 
0.000 0.000 = beta_min, beta_max ; NBETA = 1 
0.000 90.000 = theta_min, theta_max; NTHETA= 3 
0.000 0.000 = phi_min, phi_max ; NPHI = 1 
Results averaged over 3 target orientations 

and 2 incident polarizations 
Qext Qabs Qsca g(l)=<cos> Qbk Qpha 

J0=1: 1.3296E-01 3.2687E-02 1.0028E-01 2.3451E-01 5.5520E-03 4.8356E-01 
J0=2: 9.0628E-02 2.3972E-02 6.6655E-02 2.6726E-01 3.4303E-03 4.1588E-01 
mean: 1.1180E-01 2.8330E-02 8.3466E-02 2.4758E-01 4.4912E-03 4.4972E-01 
Qpol= 4.2337E-02 dQpha= 6.7685E-02 



Qsca*g(l) Qsca*g(2) Qsca*g(3) iter mxiter 



J0= 


= 1 : 
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-02 1 . 
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C ^fjxxryykzzz .sea Files 

The wOOrOOkOOO. sea file contains the results for the first wavelength (wO 0), first target radius (r 0), 
and first orientation (kO 00). The wOOrOOkOOO . sea file created by the sample calculation should look 
Uke the following: 



DDSCAT DDSCAT.6.0 [03.07.13] 

TARGET Rectangular prism; NX,NY,NZ= £ 

GPFAFT method of solution 

LATTDR prescription for polarizabilies 

RCTNGL shape 

192 = NATO = number of dipoles 



AEFF= 1.00000 = effective radius (physical units) 
WAVE= 6.28319 = wavelength (physical units) 
K*AEFF= 1.00000 = 2*pi*aeff /lambda 

n= ( 1.3300 , 0.0100), eps.= ( 1.7688 , 0.0266) |m|kd= 
TOL= l.OOOE-05 = error tolerance for CCG method 



0.3716 for subs . 



ICTHM= 33 = theta values used in comp. of Qsca, g 
IPHIM= 12 = phi values used in comp. of Qsca,g 
( 1.00000 0.00000 0.00000) = target axis Al in Target Frame 
( 0.00000 1.00000 0.00000) = target axis A2 in Target Frame 
( 0.27942 0.00000 0.00000) = k vector (latt. units) in TF 
( 0.00000, 0.00000) ( 1.00000, 0.00000) ( 0.00000, 0.00000)=inc.pol.vec. 1 
( 0.00000, 0.00000) ( 0.00000, 0.00000) ( 1.00000, 0.00000)=inc.pol.vec. 2 
BETA = 0.000 = rotation of target around Al 
THETA= 0.000 = angle between Al and k 
PHI 



in 
in 



TF 
TF 



J0=1 
J0=2 
mean 
Qpol 

J0=1 
J0=2 
mean 



0.000 

0.000 

0.000 
Qext 
1 . 1104E-01 
8 . 6508E-02 
9. 8776E-02 
2.4537E-02 

Qsca*g (1) 
2 . 8296E-02 
2 . 2663E-02 
2.5480E-02 



rotation of Al around k 



Qabs 

3.0279E-02 E 

?.4413E-02 t 

?.7345E-02 ' 

Qsca*g (2) 
-5.5908E-10 
3.0799E-09 
1 .2604E-09 



Qsca 
. 0766E-02 

, 2095E-02 
. 1430E-02 

Qsca*g (3) 
-1 .2161E-09 
1 . 4463E-09 
1.1509E-10 



g (1) =<cos> 
3.5035E-01 

3 . 6498E-01 
3 . 5671E-01 



Qbk Qpha 
1.8580E-03 4.6680E-01 

1.3620E-03 4.1968E-01 
1.6100E-03 4.4324E-01 
dQpha= 4.7115E-02 



iter 
3 
3 



mxiter 
576 
576 



** Mueller matrix elements for selected scattering directions ** 



theta 


phi 




Pol . 


c 


_11 






S_12 






S_21 






S_2 2 




c 


_31 






S_41 



















10772 


8 . 


312E- 


03 


8 


. 954E- 


04 


8 


. 954E- 


04 


8 


. 312E- 


03 


-1 . 


837E- 


11 


5 


. 826E- 


12 


30 














02296 


6. 


734E- 


03 


-1 


.546E- 


04 


-1 


.546E- 


04 


6 


. 734E- 


03 


7 . 


358E- 


13 


9 


. 139E- 


12 


60 














48383 


3. 


651E- 


03 


-1 


.766E- 


03 


-1 


.766E- 


03 


3 


. 651E- 


03 


-5. 


798E- 


12 


-1 


. 051E- 


11 


90 














99794 


1 . 


764E- 


03 


-1 


.760E- 


03 


-1 


.760E- 


03 


1 


. 764E- 


03 


-6. 


232E- 


12 


7 


. 791E- 


12 


120 














53222 


1 . 


251E- 


03 


-6 


. 658E- 


04 


-6 


. 658E- 


04 


1 


. 251E- 


03 


-1 . 


813E- 


12 


-4 


. 691E- 


13 


150 














00167 


9. 


868E- 


04 


-1 


. 644E- 


06 


-1 


. 644E- 


06 


9 


. 868E- 


04 


-5. 


853E- 


12 


-2 


.595E- 


12 


180 














15403 


8 . 


430E- 


04 


1 


.298E- 


04 


1 


.298E- 


04 


8 


.430E- 


04 


3. 


319E- 


13 


-3 


. 020E- 


12 








90 








10772 


8. 


312E- 


03 


-8 


. 954E- 


04 


-8 


. 954E- 


04 


8 


. 312E- 


03 


-7 . 


445E- 


11 


9 


. 616E- 


12 


30 





90 








24082 


7 . 


141E- 


03 


-1 


. 720E- 


03 


-1 


.720E- 


03 


7 


. 141E- 


03 


-7 . 


673E- 


11 


-6 


. 021E- 


12 


60 





90 








64798 


4 . 


560E- 


03 


-2 


. 955E- 


03 


-2 


. 955E- 


03 


4 


. 560E- 


03 


-2 . 


574E- 


11 


3 


.576E- 


12 


90 





90 








99942 


2 . 


566E- 


03 


-2 


.565E- 


03 


-2 


.565E- 


03 


2 


.566E- 


03 


-5. 


lllE- 


12 


-5 


. 125E- 


12 


120 





90 








68918 


1 . 


631E- 


03 


-1 


. 124E- 


03 


-1 


.124E- 


03 


1 


. 631E- 


03 


8. 


146E- 


12 


-8 


.819E- 


13 


150 





90 








28491 


1 . 


065E- 


03 


-3 


. 034E- 


04 


-3 


. 034E- 


04 


1 


.065E- 


03 


9. 


598E- 


12 


1 


.278E- 


12 


180 





90 








15403 


8. 


430E- 


04 


-1 


.298E- 


04 


-1 


.298E- 


04 


8 


. 430E- 


04 


7 . 


621E- 


12 


-1 


. 666E- 


12 



46 



